WSL/SLF GitLab Repository

EnergyBalance.cc 6.83 KB
Newer Older
Julien Esseiva's avatar
Julien Esseiva committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
/***********************************************************************************/
/*  Copyright 2009-2015 WSL Institute for Snow and Avalanche Research    SLF-DAVOS      */
/***********************************************************************************/
/* This file is part of Alpine3D.
    Alpine3D is free software: you can redistribute it and/or modify
    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.

    Alpine3D is distributed in the hope that it will be useful,
    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
    along with Alpine3D.  If not, see <http://www.gnu.org/licenses/>.
*/
#include <alpine3d/ebalance/EnergyBalance.h>
#include <alpine3d/MPIControl.h>
#include <alpine3d/OMPControl.h>

using namespace mio;
using namespace std;

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

Julien Esseiva's avatar
Julien Esseiva committed
32
33
34
35
36
37
38
39
40
	MPIControl& instance = MPIControl::instance();

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

	for (size_t ii=0; ii<nbworkers; ii++) {
		size_t thread_startx, thread_nx;
		OMPControl::getArraySliceParams(nx, nbworkers, ii, thread_startx, thread_nx);
		const size_t offset = startx + thread_startx;
41
42
		std::cout << "[i] EnergyBalance worker " << ii << " on process " << instance.rank() << " will start at offset " <<
		offset << " with nx " << thread_nx << "\n";
43
		radfields.push_back(RadiationField(dem_in, offset, thread_nx));
Julien Esseiva's avatar
Julien Esseiva committed
44
45
	}

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

Julien Esseiva's avatar
Julien Esseiva committed
49
50
51
	// Every MPI process will have its own copy of terrain_radiation object with full DEM
	const bool enable_terrain_radiation = cfg.get("Terrain_Radiation", "EBalance");
	if (enable_terrain_radiation) {
52
		terrain_radiation = TerrainRadiationFactory::getAlgorithm(cfg, dem, nbworkers);
Julien Esseiva's avatar
Julien Esseiva committed
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
		const std::string algo = terrain_radiation->algo;
		if (instance.master())
			std::cout << "[i] Using terrain radiation with model: " << algo << "\n";
	}
}

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

EnergyBalance& EnergyBalance::operator=(const EnergyBalance& source) {
	if (this != &source) {
		snowpack = source.snowpack;
		terrain_radiation = source.terrain_radiation;
		radfields = source.radfields;
		dem = source.dem;
		vecMeteo = source.vecMeteo;
		albedo = source.albedo;
		direct = source.direct;
		diffuse = source.diffuse;
		reflected = source.reflected;
74
		direct_unshaded_horizontal = source.reflected;
Julien Esseiva's avatar
Julien Esseiva committed
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
		timer = source.timer;
		dimx = source.dimx;
		dimy = source.dimy;
		nbworkers = source.nbworkers;
	}
	return *this;
}

std::string EnergyBalance::getGridsRequirements() const
{
	return "TOP_ALB";
}

void EnergyBalance::Destroy()
{
	if (terrain_radiation) {
		delete terrain_radiation;
		terrain_radiation = NULL;
	}
}

void EnergyBalance::setSnowPack(SnowpackInterface& mysnowpack)
{
	snowpack = &mysnowpack;
}

void EnergyBalance::setAlbedo(const mio::Grid2DObject& in_albedo)
{
	albedo = in_albedo;
104
105
106
107
108
	//resetting these grids that are not valid anymore
	direct_unshaded_horizontal=0;
	direct=0;
	diffuse=0;
	reflected=0;
Julien Esseiva's avatar
Julien Esseiva committed
109
110
111
112
113
}

void EnergyBalance::setStations(const std::vector<mio::MeteoData>& in_vecMeteo)
{
	vecMeteo = in_vecMeteo;
114
115
	//resetting these grids that are not valid anymore
	direct_unshaded_horizontal=0;
116
117
118
	direct=0;
	diffuse=0;
	reflected=0;
Julien Esseiva's avatar
Julien Esseiva committed
119
120
121
122
123
124
125
126
127
}

void EnergyBalance::setMeteo(const mio::Grid2DObject& in_ilwr,
                             const mio::Grid2DObject& in_ta, const mio::Grid2DObject& in_rh, const mio::Grid2DObject& in_p, const mio::Date timestamp)
{
	timer.restart();

	#pragma omp parallel for schedule(dynamic)
	for (size_t ii=0; ii<nbworkers; ii++) {
128
		radfields[ii].setStations(vecMeteo, albedo); //calculate the parameters at the radiation stations
Julien Esseiva's avatar
Julien Esseiva committed
129
		size_t startx, nx;
130
131
		radfields[ii].getBandOffsets(startx, nx);
		radfields[ii].setMeteo(mio::Grid2DObject(in_ta, startx, 0, nx, dimy),
Julien Esseiva's avatar
Julien Esseiva committed
132
133
134
135
		                       mio::Grid2DObject(in_rh, startx, 0, nx, dimy),
		                       mio::Grid2DObject(in_p, startx, 0, nx, dimy),
		                       mio::Grid2DObject(albedo, startx, 0, nx, dimy));

136
		mio::Array2D<double> band_direct, band_diffuse, band_direct_unshaded_horizontal;
137
		radfields[ii].getRadiation(band_direct, band_diffuse, band_direct_unshaded_horizontal);
Julien Esseiva's avatar
Julien Esseiva committed
138
139
		direct.fill(band_direct, startx, 0, nx, dimy);
		diffuse.fill(band_diffuse, startx, 0, nx, dimy);
140
		direct_unshaded_horizontal.fill(band_direct_unshaded_horizontal, startx, 0, nx, dimy);
Julien Esseiva's avatar
Julien Esseiva committed
141
142
143
	}
	MPIControl::instance().allreduce_sum(direct);
	MPIControl::instance().allreduce_sum(diffuse);
144
	MPIControl::instance().allreduce_sum(direct_unshaded_horizontal);
145
	double solarAzimuth, solarElevation;
146
	radfields[0].getPositionSun(solarAzimuth, solarElevation);
147

148
149
150
151
	if (hasSP())
		terrain_radiation->setSP(radfields[0].getDate(), solarAzimuth, solarElevation);


Julien Esseiva's avatar
Julien Esseiva committed
152

153
	mio::Array2D<double> view_factor(dimx, dimy, IOUtils::nodata);
Julien Esseiva's avatar
Julien Esseiva committed
154
155
156
	if (terrain_radiation) {
		// note: parallelization has to take place inside the TerrainRadiationAlgorithm implementations
		terrain_radiation->setMeteo(albedo.grid2D, in_ta.grid2D, in_rh.grid2D, in_ilwr.grid2D);
157
158
		terrain_radiation->getRadiation(direct, diffuse, reflected, direct_unshaded_horizontal,view_factor,
                                    solarAzimuth, solarElevation);
Julien Esseiva's avatar
Julien Esseiva committed
159
160
161
162
163
164
165
166
167
168
169
170
171
	}

	if (MPIControl::instance().master())
		cout << "[i] Ebalance simulation done for " << timestamp.toString(Date::ISO) << "\n";

	if (snowpack) {
		mio::Array2D<double> ilwr = in_ilwr.grid2D;
		mio::Array2D<double> global = direct+diffuse; //otherwise the compiler does not match the types

		if (!reflected.empty()) global += reflected;

		timer.stop();
		try {
172
173
			snowpack->setRadiationComponents(global, ilwr, diffuse, view_factor, reflected,
                                       in_ilwr.grid2D,  solarElevation, timestamp); //this triggers Snowpack calculation
Julien Esseiva's avatar
Julien Esseiva committed
174
175
176
177
178
179
180
181
182
		} catch(std::exception& e) {
			std::cout << "[E] Exception in snowpack->setRadiationComponents()\n";
			cout << e.what() << endl;
			std::abort(); //force core dump
		}
	}
	timer.stop();
}

Adrien Michel's avatar
Adrien Michel committed
183
void EnergyBalance::writeSP(const unsigned int max_steps){
184
	if (hasSP()) terrain_radiation->writeSP(max_steps);
185
186
}

187

Julien Esseiva's avatar
Julien Esseiva committed
188
189
190
191
double EnergyBalance::getTiming() const
{
	return timer.getElapsed();
}