Commit 6f504fd2 authored by Nicolas Rodriguez's avatar Nicolas Rodriguez
Browse files

Upload New File

parent 677bc36f
% Copyright 2019-2020 LIST (Luxembourg Institute of Science and Technology), all right reserved.
% Author: Nicolas Rodriguez (
% This is the main script calling all the other scripts in order:
% init_virtual initializes the model: parameters, constants, and loads hydrological and tracer data
% model_run_virtual governs the model runs, calling the calculation script model_calc, and evaluating model performance
% The last part of the script is the Monte Carlo for parameter calibration
% Customizations should be mainly specified in this script
% Clear and close
clear variables;
close all;
% Start clock
% Display options (text display about calculation status in the command window)
options.display=1; % 0 for no display, 1 for all display
options.plots.display=0; % 0 for no plots, 1 for all plots
% Solver options (numerical solution to water mass balance equation)
options.solver.type='EE'; % Explicit Euler, 'IE' for implicit Euler
options.solver.checks=0.25; % How often to check calculations? (as a percentage of completion)
if strcmp(options.solver.type,'IE') % Implicit Euler, which is more or less useless in this context
% see file mysolver.m for the parameters below
options.solver.tol=eps; % approximation of zero for function evaluations
options.solver.prec=eps; % convergence criteria on the difference between function evaluations
options.solver.maxiter=1000; % max number of iterations for the solver
options.solver.display='on'; % 'off' for no solver messages
options.solver.maxdisp=100; % max nb of msgs by solver
% Model parameter options
% SASQ_p1 to SASQ_pn: n parameters for the streamflow SAS function
% same for ET SAS function
% other parameters (do not need to be changed):
% m: smoothing parameter for AET calculations from PET (see supplement info and model_calc.m)
% S_ET: storage value below which AET drops from PET towards 0, using the smoothing parameter m
% n: smoothing parameter for storage calculated from WB (see model_calc.m)
% Psi: optional multiplicative constant for ET
% C_i: initial tracer concentration in storage (per mil)
% S_i: initial storage (mm)
% T0: mean age in years for the initial residence time distribution (exponential)
% Calibration parameters (if needed, more parameters can be taken from the fixed parameter list and put here, for instance, one can cut 'SASQ_p3' and paste here, and adjust options.param.fix.nb to 22 and to 3; % number of calibration parameters (here only two are needed if the unimodal gamma SAS function is used, they will correspond to SASQ_p1 = k and SASQ_p2 = theta){'SASQ_p1';'SASQ_p2'};
% p calibration parameters are placed in a vector param(1,p)
% Fixed parameters (their values can be modified only in init_virtual.m)
options.param.fix.nb=23; % number of fixed parameters (not calibrated){'m';'T0';'n';'S_i';'C_i';'Psi';'S_ET';'SASQ_p3';'SASQ_p4';'SASQ_p5';'SASQ_p6';'SASQ_p7';'SASQ_p8';'SASQ_p9';'SASET_p1';'SASET_p2';'SASET_p3';'SASET_p4';'SASET_p5';'SASET_p6';'SASET_p7';'SASET_p8';'SASET_p9'};
% all other parameters are placed in the 'options' stucture, and fixed during initialisation
% Age calculation options
options.ages.dT=24*2; % age resolution for the calculations, in hours
options.ages.lookup=1; % use a table lookup for the gamma distributions (instead of the built-in function gampdf)? This can make calculations faster
options.ages.dS=1; % resolution of age-ranked storage used to lookup the gamma pdfs
options.ages.short=1; % shorten the TTDs by neglecting numbers smaller than eps in the tails? This can make calculations faster
options.ages.norm=0; % normalize distributions to avoid numerical inaccruacies? (basically make sure the weights applied in the convolution sum to 1)
options.ages.save_warmup=0; % save age results of the warmup? (uses a lot of memory!)
options.ages.save_interval=24*40; % interval of time (in the simulation time) to save age calculations, in hours (needed to reduce memory and RAM constraints)
options.ages.init_RTD=[]; % initial residence time distribution in the reservoir that can be loaded instead of the default exponential distribution with mean T0
options.ages.age_prctiles=[0.10 0.25 0.50 0.75 0.90]; % percentiles of TTDs and RTDs to be calculated
options.ages.T_group=24*[2 7 31 50 73 100 120 150 240 360 365 550 1000 1100 3000]; % groups to classify water younger than the given ages hours (needs to be round number)
options.ages.SAStype_Q='gam'; % type of SAS function used for Q, see evalSAS.m (here gamma distribution)
options.ages.SAStype_ET='RS'; % type of SAS function used for ET, see evalSAS.m (here Random Sampling)
% Tracer options
options.tracer=1; % 1 for 2H (only available choice here)
% Initialization scripts
options.dt=24; % desired data time step in hours
options.y0=[]; % no user-loaded initial values for state variables
options.nbyrs_warmup=0; % number of years for the warmup period (here there is no warmup needed, we directly match the simulated values to observations that already exclude the 12.5 first years as warmup)
% Run the model
options.movavg=1; % apply moving average to results? (necessary to remove artificial numerical oscillations due to the convolution equation and the different resolutions dt and dT);
param=[624/87,87]; % k=mu/theta and theta of the SAS of Q gamma distribution (best set from the monte carlo)
[Eff_C,model_output,model_output_warmup,options]=model_run_virtual(param, model_input, meas, options); % script that runs the model
% model_output contains all simulation results, including all age results (see model_calc.m)
%% Monte Carlo script
% Note: this was run on a high performance computer (using a parfor loop). A standard office PC would take about 15 days too complete the Monte Carlo simulation!
% The results can be found in the file "Monte_Carlo_results.mat". Columns 1 and 2 respectively correspond to parameters k=mu/theta and theta, column 3 corresponds to the NSE
% sz=15000; % number of parameter sets
% results_monte_carlo=zeros(sz,3);
% for k=1:sz
% mu=unifrnd(1,700);
% theta=unifrnd(1,500);
% param=[mu/theta, theta];
% [Eff_C,~,~,~]=model_run_virtual(param, model_input, meas, options);
% NS=Eff_C;
% results_monte_carlo(k,:)=cat(2,param,NS);
% end
disp(char(['Total time spent: ' num2str(t2) ' s']));
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment