WSL/SLF GitLab Repository

EnergyBalance.cc 6.23 KB
Newer Older
1
/***********************************************************************************/
2
/*  Copyright 2009-2015 WSL Institute for Snow and Avalanche Research    SLF-DAVOS      */
3
/***********************************************************************************/
4
5
/* This file is part of Alpine3D.
    Alpine3D is free software: you can redistribute it and/or modify
6
7
8
9
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

10
    Alpine3D is distributed in the hope that it will be useful,
11
12
13
14
15
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    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
#include <alpine3d/ebalance/EnergyBalance.h>
19

20
21
22
using namespace mio;
using namespace std;

23
EnergyBalance::EnergyBalance(const unsigned int& i_nbworkers, const mio::Config& cfg, const mio::DEMObject &dem_in)
24
25
26
              : snowpack(NULL), terrain_radiation(NULL), radfields(i_nbworkers), dem(dem_in), vecMeteo(),
                albedo(dem_in, 0.), direct(), diffuse(), reflected(),
                timer(), dimx(dem_in.getNx()), dimy(dem_in.getNy()), nbworkers(i_nbworkers)
27
{
28
29
30
31
32
	MPIControl& instance = MPIControl::instance();

	size_t startx = 0, nx = dimx;
	instance.getArraySliceParams(dimx, startx, nx);

33
	#pragma omp parallel for schedule(static)
34
35
36
	for (size_t ii=0; ii<nbworkers; ii++) {
		size_t thread_startx, thread_nx;
		IOUtils::getArraySliceParams(nx, nbworkers, ii+1, thread_startx, thread_nx);
37
		const size_t offset = startx + thread_startx;
38

39
		#pragma omp critical(ebWorkers_status)
40
		std::cout << "[i] EnergyBalance worker " << ii << " on process " << instance.rank() << " will start at offset " << offset << " with nx " << thread_nx << "\n";
41
		radfields[ii] = new RadiationField(dem_in, offset, thread_nx);
42
	}
43

44
	if (instance.master())
45
		std::cout << "[i] EnergyBalance initialized a total of " << instance.size() << " process(es) with " << nbworkers << " worker(s) each\n";
46

47
	// Every MPI process will have its own copy of terrain_radiation object with full DEM
48
	const bool enable_terrain_radiation = cfg.get("Terrain_Radiation", "EBalance");
49
	if (enable_terrain_radiation) {
50
		terrain_radiation = TerrainRadiationFactory::getAlgorithm(cfg, dem, nbworkers);
51
		const std::string algo = terrain_radiation->algo;
52
53
		if (instance.master())
			std::cout << "[i] Using terrain radiation with model: " << algo << "\n";
54
	}
55
56
57
58
59
60
}

EnergyBalance::~EnergyBalance() {
	Destroy( );
}

61
62
63
64
65
66
EnergyBalance& EnergyBalance::operator=(const EnergyBalance& source) {
	if (this != &source) {
		snowpack = source.snowpack;
		terrain_radiation = source.terrain_radiation;
		radfields = source.radfields;
		dem = source.dem;
67
68
		vecMeteo = source.vecMeteo;
		albedo = source.albedo;
69
70
71
72
73
74
75
76
77
78
79
		direct = source.direct;
		diffuse = source.diffuse;
		reflected = source.reflected;
		timer = source.timer;
		dimx = source.dimx;
		dimy = source.dimy;
		nbworkers = source.nbworkers;
	}
	return *this;
}

80
81
std::string EnergyBalance::getGridsRequirements() const
{
82
	return "TOP_ALB";
83
84
}

85
86
void EnergyBalance::Destroy()
{
87
88
89
90
91
92
93
94
95
	while (!radfields.empty()) {
		delete radfields.back();
		radfields.pop_back();
	}

	if (terrain_radiation) {
		delete terrain_radiation;
		terrain_radiation = NULL;
	}
96
97
}

Thomas Egger's avatar
Thomas Egger committed
98
void EnergyBalance::setSnowPack(SnowpackInterface& mysnowpack)
99
100
101
102
{
	snowpack = &mysnowpack;
}

Thomas Egger's avatar
Thomas Egger committed
103
104
void EnergyBalance::setAlbedo(const mio::Grid2DObject& in_albedo)
{
105
	albedo = in_albedo;
106
107
108

	direct.resize(0, 0); //resetting these grids that are not valid anymore
	diffuse.resize(0, 0);
109
	reflected.resize(0, 0);
110
111
}

112
void EnergyBalance::setStations(const std::vector<mio::MeteoData>& in_vecMeteo)
Thomas Egger's avatar
Thomas Egger committed
113
{
114
	vecMeteo = in_vecMeteo;
115
116
117

	direct.resize(0, 0); //resetting these grids that are not valid anymore
	diffuse.resize(0, 0);
118
	reflected.resize(0, 0);
119
120
121
}

void EnergyBalance::setMeteo(const mio::Grid2DObject& in_ilwr,
Thomas Egger's avatar
Thomas Egger committed
122
123
                             const mio::Grid2DObject& in_ta, const mio::Grid2DObject& in_rh, const mio::Grid2DObject& in_p, const mio::Date timestamp)
{
124
	timer.restart();
125
126
	direct.resize(dimx, dimy);
	diffuse.resize(dimx, dimy);
127
	
128
	#pragma omp parallel for schedule(dynamic)
129
	for (size_t ii=0; ii<nbworkers; ii++) {
130
		radfields[ii]->setStations(vecMeteo, albedo); //calculate the parameters at the radiation stations
131
132
		size_t startx, nx;
		radfields[ii]->getBandOffsets(startx, nx);
133
134
135
136
		radfields[ii]->setMeteo(mio::Grid2DObject(in_ta, startx, 0, nx, dimy),
		                       mio::Grid2DObject(in_rh, startx, 0, nx, dimy),
		                       mio::Grid2DObject(in_p, startx, 0, nx, dimy),
		                       mio::Grid2DObject(albedo, startx, 0, nx, dimy));
137
138
		
		mio::Array2D<double> band_direct, band_diffuse;
139
140
141
142
143
144
		radfields[ii]->getRadiation(band_direct, band_diffuse);
		direct.fill(band_direct, startx, 0, nx, dimy);
		diffuse.fill(band_diffuse, startx, 0, nx, dimy);
	}
	MPIControl::instance().allreduce_sum(direct);
	MPIControl::instance().allreduce_sum(diffuse);
145
146

	if (terrain_radiation) {
147
		// note: parallelization has to take place inside the TerrainRadiationAlgorithm implementations
148
		terrain_radiation->setMeteo(albedo.grid2D, in_ta.grid2D, in_rh.grid2D, in_ilwr.grid2D);
149
		terrain_radiation->getRadiation(direct, diffuse, reflected);
150
	}
151

Thomas Egger's avatar
Thomas Egger committed
152
	if (MPIControl::instance().master())
153
		cout << "[i] Ebalance simulation done for " << timestamp.toString(Date::ISO) << "\n";
154

155
	if (snowpack) {
156
157
158
		double solarAzimuth, solarElevation;
		radfields[0]->getPositionSun(solarAzimuth, solarElevation); //we need it only for handing over to snowpack
		
159
160
		mio::Array2D<double> ilwr = in_ilwr.grid2D;
		mio::Array2D<double> global = direct+diffuse; //otherwise the compiler does not match the types
161

162
		if (!reflected.empty()) global += reflected;
163
		
164
		timer.stop();
165
166
167
		try {
			snowpack->setRadiationComponents(global, ilwr, diffuse, solarElevation, timestamp); //this triggers Snowpack calculation
		} catch(std::exception& e) {
168
			std::cout << "[E] Exception in snowpack->setRadiationComponents()\n";
169
			cout << e.what() << endl;
170
			std::abort(); //force core dump
171
		}
172
173
174
175
176
177
178
179
180
	}
	timer.stop();
}

double EnergyBalance::getTiming() const
{
	return timer.getElapsed();
}