Alpine3D

Alpine3D Svn Source Tree

Root/trunk/alpine3d/AlpineMain.cc

1/***********************************************************************************/
2/* Copyright 2009-2015 WSL Institute for Snow and Avalanche Research SLF-DAVOS */
3/***********************************************************************************/
4/* This file is part of Alpine3D.
5 Alpine3D is free software: you can redistribute it and/or modify
6 it under the terms of the GNU Lesser General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9
10 Alpine3D is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU Lesser General Public License for more details.
14
15 You should have received a copy of the GNU Lesser General Public License
16 along with Alpine3D. If not, see <http://www.gnu.org/licenses/>.
17*/
18#ifdef DEBUG_ARITHM
19#ifndef _GNU_SOURCE
20#define _GNU_SOURCE
21#endif
22#include <fenv.h>
23#endif
24
25#include <exception>
26#include <iostream>
27#include <csignal>
28#include <unistd.h>
29#include <getopt.h>
30#include <ctime>
31
32#include <meteoio/MeteoIO.h>
33#include <snowpack/libsnowpack.h>
34
35#include <alpine3d/AlpineMain.h>
36#include <alpine3d/AlpineControl.h>
37#include <alpine3d/MPIControl.h>
38
39#include <alpine3d/ebalance/EnergyBalance.h>
40#include <alpine3d/snowdrift/SnowDrift.h>
41#include <alpine3d/SnowpackInterface.h>
42#include <alpine3d/DataAssimilation.h>
43
44using namespace std;
45using namespace mio;
46
47//Global variables for AlpineMain
48static int steps = 0; //how many hours to simulate
49static Date startdate;
50
51static bool enable_eb=false, enable_drift=false, enable_runoff=false, enable_da=false, nocompute=false, restart=false;
52static int npsnowpack=1, npebalance=1;
53
54inline void Usage(char *prog)
55{
56cout << "Usage: " << prog << "\n";
57cout << "\t-a, --startdate=YYYY-MM-DDTHH:MM (e.g.:2008-08-11T09:00)\n";
58cout << "\t-z, --enddate=YYYY-MM-DDTHH:MM (e.g.:2008-08-11T09:00) OR ";
59cout << "-n, --steps=<number of timesteps>\n";
60cout << "\t[--enable-eb] (default:off)\n";
61cout << "\t[--enable-runoff] (default:off)\n";
62cout << "\t[--enable-drift] (default:off)\n";
63cout << "\t[--enable-da] (default:off)\n";
64cout << "\t[--no-compute] Don't compute the simulation, only parse and validate the inputs (default:off)\n";
65cout << "\t[--restart] Restart simulation from previous .sno files (default:off)\n";
66cout << "\t[-p, --np-snowpack=<number of processors for SnowPack>]\n";
67cout << "\t[-b, --np-ebalance=<number of processors for Ebalance>]\n";
68cout << "\t[-i, --iofile=<MeteoIO ini file> (default:./io.ini)]\n";
69cout << "\t[-h, --help] Print help message and version information\n";
70
71exit(1);
72}
73
74inline void parseCmdLine(int argc, char **argv, Config &cfg)
75{
76std::string iofile( "io.ini" );
77int eeb=0, edr=0, ero=0, eda=0, nco=0, rst=0;
78int longindex=0, opt = -1;
79bool setStart=false, setEnd=false, setSteps=false;
80std::string start_date_input, end_date_input;
81
82const struct option long_options[] =
83{
84/* These options set a flag. */
85{"enable-eb",no_argument,&eeb, 1},
86{"enable-drift",no_argument,&edr, 1},
87{"enable-runoff",no_argument,&ero, 1},
88{"enable-da",no_argument,&eda, 1},
89{"no-compute",no_argument,&nco, 1},
90{"restart",no_argument,&rst, 1},
91
92/* These options don't set a flag.
93We distinguish them by their indices. */
94{"np-snowpack",required_argument, 0, 'p'},
95{"np-ebalance",required_argument, 0, 'b'},
96{"iofile",required_argument, 0, 'i'},
97{"startdate",required_argument, 0, 'a'},
98{"enddate",required_argument, 0, 'z'},
99{"steps",required_argument, 0, 'n'},
100{"help",no_argument, 0, 'h'},
101{0, 0, 0, 0}
102};
103
104while ((opt=getopt_long( argc, argv, ":d:p:r:s:e:i:f:a:z:n:h", long_options, &longindex)) != -1) {
105switch (opt){
106case 0:
107break;
108case 'b':
109IOUtils::convertString(npebalance, string(optarg));
110if (npebalance < 1)
111throw IOException("Number of processors for Ebalance computation is not valid!", AT);
112break;
113case 'p':
114IOUtils::convertString(npsnowpack, string(optarg));
115if (npsnowpack < 1)
116throw IOException("Number of processors for SnowPack computation is not valid!", AT);
117break;
118case 'i':
119iofile = string(optarg);
120break;
121case 'a':
122start_date_input = string(optarg);
123setStart = true;
124break;
125case 'z':
126end_date_input = string(optarg);
127setEnd = true;
128break;
129case 'n':
130IOUtils::convertString(steps, string(optarg));
131setSteps = true;
132break;
133
134case ':': //operand missing
135cout << "Command line parameter '-" << char(optopt) << "' requires an operand";
136break;
137case 'h':
138Usage(argv[0]);
139break;
140case '?':
141cerr << "\nERROR: Unknown argument detected\n";
142Usage(argv[0]);
143break;
144default:
145cerr << "\nERROR: getopt returned character code " << opt << "\n";
146Usage(argv[0]);
147}
148}
149
150//check that the minimum flags have been provided by the user
151if (!((setStart && setSteps) || (setStart && setEnd)) && MPIControl::instance().master()) {
152cout << "\nERROR: You must at least specify the 'startdate' and the 'steps' parameters"
153 << " or the 'startdate' and the 'enddate' parameters\n\n";
154Usage(argv[0]);
155}
156
157enable_eb= (eeb!=0);
158enable_drift= (edr!=0);
159enable_runoff= (ero!=0);
160enable_da= (eda!=0);
161nocompute= (nco!=0);
162restart= (rst!=0);
163
164//now, we read the config file and set the start and end dates
165//only the master requires the io.ini, the rest receives a the cfg object:
166const bool isMaster = MPIControl::instance().master(); // Check if this process is the master (always true for non-parallel mode)
167if (isMaster) {
168cfg.addFile(iofile);
169}
170MPIControl::instance().broadcast(cfg);
171
172const double TZ = cfg.get("TIME_ZONE", "INPUT");
173startdate.setTimeZone(TZ);
174Date enddate;
175enddate.setTimeZone(TZ);
176IOUtils::convertString(startdate, start_date_input, TZ);
177if (end_date_input == "NOW") { //interpret user provided end date
178enddate.setFromSys();
179enddate.setTimeZone(TZ);
180enddate.rnd(1800, Date::DOWN);
181} else {
182IOUtils::convertString(enddate, end_date_input, TZ);
183}
184
185//check that the start and end dates make sense
186if (setStart && setEnd) {
187steps = (int)floor((enddate.getJulian(true) - startdate.getJulian(true))*24.+0.01) + 1; //to be on the safe side when rounding
188if (steps <= 0) {
189 throw InvalidArgumentException("startdate="+startdate.toString(Date::ISO)+" > enddate="+enddate.toString(Date::ISO),AT);
190}
191}
192
193//SNOWPARAM: Check for time step consistency
194const double sp_dt = cfg.get("CALCULATION_STEP_LENGTH", "Snowpack");
195if ( dt_main != M_TO_S(sp_dt) ) {
196if (isMaster) cout << "[I] Alpine3D's dt_main=" << S_TO_M(dt_main) << " [min], Snowpack's CALCULATION_STEP_LENGTH=" << sp_dt << " [min] (in io.ini)\t";
197
198if ( (int)(dt_main)%(int)(M_TO_S(sp_dt)) == 0) {
199if (isMaster) cout << "--> " << (int)(dt_main/M_TO_S(sp_dt)) << " Snowpack simulations per dt_main\n";
200} else {
201if (isMaster) cout << "Time step inconsistency between Alpine3D's dt_main=" << S_TO_M(dt_main) << " [min], Snowpack's CALCULATION_STEP_LENGTH=" << sp_dt << " [min] (in io.ini). Please set Snowpack's time step to an integral fraction of dt_main!\n";
202exit(1);
203}
204}
205}
206
207inline void setStaticData(const Config &cfg, IOManager& io, DEMObject &dem, Grid2DObject &landuse, std::vector<Coords> &vec_pts)
208{
209const bool isMaster = MPIControl::instance().master(); // Check if this process is the master (always true for non-parallel mode)
210
211if (isMaster) { //Reading DEM, LANDUSE, special PTS on master process only
212const std::string slope_algorithm = cfg.get("SLOPE_ALGORITHM", "Input", "CORRIPIO");
213dem.setDefaultAlgorithm(slope_algorithm);
214dem.setUpdatePpt((DEMObject::update_type)(DEMObject::SLOPE | DEMObject::NORMAL | DEMObject::CURVATURE));
215
216io.readDEM(dem);
217
218//cleanup points that might have elevation but no slope: this makes for much easier checks!
219dem.sanitize();
220io.readLanduse(landuse);
221//TODO: generate a map of bool saying for each pixel if it is a valid pixel or not
222//-> if (isValid(ix,iy)==true) -> calculate (eb, sn, etc)
223cerr << "[i] Read 2D Grid DEM: " << dem.getNx() << "x" << dem.getNy() << "\n";
224cerr << "[i] Read 2D Grid Landuse: " << landuse.getNx() << "x" << landuse.getNy() << "\n";
225if (!landuse.isSameGeolocalization(dem))
226throw IOException("The landuse and the provided DEM don't have the same geolocalization!", AT);
227
228if (cfg.keyExists("POI", "Input")) {
229io.readPOI(vec_pts);
230if (!dem.gridify(vec_pts, true)) { //keep invalid points
231if (isMaster) cerr << "[E] Some POI are invalid or outside the DEM:\n";
232for (size_t ii=0; ii<vec_pts.size(); ii++)
233if (!vec_pts[ii].indexIsValid() && isMaster)
234std::cout << "[E] Point " << ii << "\t" << vec_pts[ii].toString(Coords::CARTESIAN) << "\n";
235throw InvalidArgumentException("Invalid POI, please check in the logs", AT);
236} else if (isMaster)
237std::cout << "[i] Using " << vec_pts.size() << " POI\n";
238}
239
240bool local_coords = false;
241cfg.getValue("COMPUTE_IN_LOCAL_COORDS", "Input", local_coords, IOUtils::nothrow);
242if (local_coords) {
243dem.llcorner.setProj("LOCAL","");
244dem.llcorner.setLocalRef(dem.llcorner.getLat(), dem.llcorner.getLon());
245}
246landuse.llcorner.copyProj(dem.llcorner);
247}
248
249MPIControl::instance().broadcast(dem);
250MPIControl::instance().broadcast(landuse);
251MPIControl::instance().broadcast(vec_pts);
252}
253
254inline void setModules(const Config &cfg, IOManager& io, const DEMObject &dem, const Grid2DObject &landuse, const std::vector<Coords> &vec_pts, SnowDriftA3D*& drift, EnergyBalance*& eb, SnowpackInterface*& snowpack, DataAssimilation*& da, Runoff*& runoff)
255{
256const bool isMaster = MPIControl::instance().master(); // Check if this process is the master (always true for non-parallel mode)
257
258//EBALANCE
259if (enable_eb && !nocompute) {
260try {
261eb = new EnergyBalance(npebalance, cfg, dem);
262} catch(std::exception& e) {
263std::cout << "[E] Exception in EnergyBalance constructor\n";
264cout << e.what() << endl;
265throw;
266}
267}
268
269//SNOWDRIFT
270if (enable_drift && !nocompute && isMaster) {
271try {
272drift = new SnowDriftA3D(dem, cfg);
273} catch(std::exception& e) {
274std::cout << "[E] Exception in SnowDriftA3D constructor\n";
275cout << e.what() << endl;
276throw;
277}
278}
279
280//DATA ASSIMILATION
281if (enable_da && !nocompute && isMaster) {
282da = new DataAssimilation(io);
283}
284
285//RUNOFF
286if (enable_runoff) {
287SnowpackConfig sn_cfg( cfg ); //so we also get the default value of "THRESH_RAIN"
288const double thresh_rain = sn_cfg.get("THRESH_RAIN", "SnowpackAdvanced");
289runoff = new Runoff (cfg, dem, IOUtils::C_TO_K(thresh_rain));
290}
291
292//SNOWPACK
293const bool glacier_katabatic_flow = cfg.get("GLACIER_KATABATIC_FLOW", "Alpine3D", false);
294if (!nocompute || glacier_katabatic_flow){
295try {
296std::string grids_requirements;
297if (enable_eb) grids_requirements = grids_requirements+" "+eb->getGridsRequirements();
298if (enable_drift) grids_requirements = grids_requirements+" "+drift->getGridsRequirements();
299if (enable_da) grids_requirements = grids_requirements+" "+da->getGridsRequirements();
300if (enable_runoff) grids_requirements = grids_requirements+" "+runoff->getGridsRequirements();
301snowpack = new SnowpackInterface(cfg, npsnowpack, dem, landuse, vec_pts, startdate, grids_requirements, restart);
302} catch(std::exception& e) {
303std::cout << "[E] Exception in SnowpackInterface constructor\n";
304cout << e.what() << endl;
305throw;
306}
307}
308
309//set callback pointers
310if (!nocompute) {
311if (drift) {
312snowpack->setSnowDrift(*drift);
313drift->setSnowPack(*snowpack);
314if (eb) drift->setEnergyBalance(*eb);
315if (isMaster) cout << "[i] Snowpack and Snowdrift interfaces exchanged\n";
316}
317
318if (eb) {
319snowpack->setEnergyBalance(*eb);
320eb->setSnowPack(*snowpack);
321if (isMaster) cout << "[i] Snowpack and Ebalance interfaces exchanged\n";
322}
323
324if (da) {
325snowpack->setDataAssimilation(*da);
326da->setSnowPack(*snowpack);
327if (isMaster) cout << "[i] Snowpack and Ebalance interfaces exchanged\n";
328}
329
330if (runoff) {
331snowpack->setRunoff(*runoff);
332runoff->setSnowPack(*snowpack);
333if (isMaster) cout << "[i] Snowpack and Ebalance interfaces exchanged\n";
334}
335}
336}
337
338inline void cleanDestroyAll(SnowDriftA3D*& drift, EnergyBalance*& eb, SnowpackInterface*& snowpack, DataAssimilation*& da, Runoff*& runoff)
339{
340if (da) da->Destroy();
341if (eb) eb->Destroy();
342if (drift) drift->Destroy();
343
344if (da) delete da;
345if (eb) delete eb;
346if (drift) delete drift;
347if (runoff) delete runoff;
348if (snowpack) delete snowpack;
349}
350
351inline void start_message(int argc, char **argv)
352{
353MPIControl& mpicontrol = MPIControl::instance();
354
355cout << argv[0] << " " << A3D_VERSION << " compiled on " << __DATE__ << " " << __TIME__ << "\n";
356cout << "\tLibsnowpack " << snowpack::getLibVersion() << "\n";
357cout << "\tMeteoIO " << mio::getLibVersion() << "\n";
358if (argc==1) Usage(argv[0]);
359
360const size_t n_process = mpicontrol.size();
361const size_t n_threads = mpicontrol.max_threads();
362cout << "\nRunning as " << n_process << " process";
363if (n_process>1) cout << "es";
364cout << " and " << n_threads << " thread";
365if (n_threads>1) cout << "s";
366cout << " with command line:";
367for (int args=0; args<argc; args++)
368cout << " " << argv[args];
369cout << endl;
370
371}
372
373inline void real_main(int argc, char **argv)
374{
375#ifdef DEBUG_ARITHM
376feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW );
377#endif
378MPIControl& mpicontrol = MPIControl::instance();
379mpicontrol.barrier(); //make sure all nodes are initialized
380if (mpicontrol.master()) {
381start_message(argc, argv);
382}
383
384Config cfg;
385parseCmdLine(argc, argv, cfg); //pass command line arguments, set config file for cfg
386IOManager io(cfg);
387const time_t start = time(NULL);
388
389//Try doing a struct "modules" with pointers on modules -> easier to pass?
390SnowDriftA3D *drift(NULL);
391EnergyBalance *eb(NULL);
392SnowpackInterface *snowpack(NULL);
393DataAssimilation* da(NULL);
394Runoff* runoff(NULL);
395
396DEMObject dem;
397Grid2DObject landuse;
398std::vector<Coords> vec_pts;
399
400//make sure that if the sno files are not written out, they will still be available for the POI
401const std::string meteo_outpath = cfg.get("METEOPATH", "Output");
402const bool snow_write = cfg.get("SNOW_WRITE", "Output");
403if (!snow_write) {
404cfg.addKey("SNOWPATH", "Output", meteo_outpath);
405}
406cfg.write(meteo_outpath + "/io.ini"); //backup the ini file
407
408try { //main integration loop
409setStaticData(cfg, io, dem, landuse, vec_pts);
410setModules(cfg, io, dem, landuse, vec_pts, drift, eb, snowpack, da, runoff);
411AlpineControl control(snowpack, drift, eb, da, runoff, cfg, dem);
412control.setNoCompute(nocompute);
413control.Run(startdate, steps);
414} catch (std::exception& e) {
415cerr << "[E] exception thrown: " << e.what() << endl;
416cleanDestroyAll(drift, eb, snowpack, da, runoff);
417exit(1);
418} catch (...) {
419cout << "[E] exception thrown!" << endl;
420cleanDestroyAll(drift, eb, snowpack, da, runoff);
421exit(1);
422}
423
424cleanDestroyAll(drift, eb, snowpack, da, runoff);
425const time_t end = time(NULL);
426if (mpicontrol.master()) {
427printf("\n");
428printf(" STARTED Alpine3D Model on %s", ctime(&start));
429printf(" ===========================================================================\n");
430printf(" FINISHED Alpine3D Model on %s", ctime(&end) );
431printf(" ===========================================================================\n");
432}
433}
434
435int main(int argc, char *argv[]) {
436try {
437real_main(argc, argv);
438} catch (const std::exception &e) {
439std::cerr << e.what() << endl;
440exit(1);
441}
442
443return 0;
444}

Archive Download this file

Revision: HEAD