WSL/SLF GitLab Repository

EnergyBalance.cc 7.2 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
		const std::string algo = terrain_radiation->algo;
		if (instance.master())
			std::cout << "[i] Using terrain radiation with model: " << algo << "\n";
	}
57
58
59
60
61
62
63
64
65
66
	bool write_sky_vf=false;
	cfg.getValue("WRITE_SKY_VIEW_FACTOR", "output", write_sky_vf,IOUtils::nothrow);

	if(MPIControl::instance().master() && write_sky_vf){
		std::cout << "[i] Writing sky view factor grid" << std::endl;
		mio::IOManager io(cfg);
		mio::Array2D<double> sky_vf(dimx,dimy,mio::IOUtils::nodata);
		terrain_radiation->getSkyViewFactor(sky_vf);
		io.write2DGrid(mio::Grid2DObject(dem_in.cellsize,dem_in.llcorner,sky_vf), "SKY_VIEW_FACTOR");
	}
Julien Esseiva's avatar
Julien Esseiva committed
67
68
}

69

Julien Esseiva's avatar
Julien Esseiva committed
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
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;
85
		direct_unshaded_horizontal = source.reflected;
Julien Esseiva's avatar
Julien Esseiva committed
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
		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;
115
116
117
118
119
	//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
120
121
122
123
124
}

void EnergyBalance::setStations(const std::vector<mio::MeteoData>& in_vecMeteo)
{
	vecMeteo = in_vecMeteo;
125
126
	//resetting these grids that are not valid anymore
	direct_unshaded_horizontal=0;
127
128
129
	direct=0;
	diffuse=0;
	reflected=0;
Julien Esseiva's avatar
Julien Esseiva committed
130
131
132
133
134
135
136
137
138
}

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++) {
139
		radfields[ii].setStations(vecMeteo, albedo); //calculate the parameters at the radiation stations
Julien Esseiva's avatar
Julien Esseiva committed
140
		size_t startx, nx;
141
142
		radfields[ii].getBandOffsets(startx, nx);
		radfields[ii].setMeteo(mio::Grid2DObject(in_ta, startx, 0, nx, dimy),
Julien Esseiva's avatar
Julien Esseiva committed
143
144
145
146
		                       mio::Grid2DObject(in_rh, startx, 0, nx, dimy),
		                       mio::Grid2DObject(in_p, startx, 0, nx, dimy),
		                       mio::Grid2DObject(albedo, startx, 0, nx, dimy));

147
		mio::Array2D<double> band_direct, band_diffuse, band_direct_unshaded_horizontal;
148
		radfields[ii].getRadiation(band_direct, band_diffuse, band_direct_unshaded_horizontal);
Julien Esseiva's avatar
Julien Esseiva committed
149
150
		direct.fill(band_direct, startx, 0, nx, dimy);
		diffuse.fill(band_diffuse, startx, 0, nx, dimy);
151
		direct_unshaded_horizontal.fill(band_direct_unshaded_horizontal, startx, 0, nx, dimy);
Julien Esseiva's avatar
Julien Esseiva committed
152
153
154
	}
	MPIControl::instance().allreduce_sum(direct);
	MPIControl::instance().allreduce_sum(diffuse);
155
	MPIControl::instance().allreduce_sum(direct_unshaded_horizontal);
156
	double solarAzimuth, solarElevation;
157
	radfields[0].getPositionSun(solarAzimuth, solarElevation);
158

159
160
161
	if (hasSP())
		terrain_radiation->setSP(radfields[0].getDate(), solarAzimuth, solarElevation);

Julien Esseiva's avatar
Julien Esseiva committed
162
163
164
	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);
165
		terrain_radiation->getRadiation(direct, diffuse, reflected, direct_unshaded_horizontal,
166
                                    solarAzimuth, solarElevation);
Julien Esseiva's avatar
Julien Esseiva committed
167
168
169
170
171
172
173
174
175
176
177
178
179
	}

	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 {
180
			snowpack->setRadiationComponents(global, ilwr, diffuse, reflected,
181
                                       in_ilwr.grid2D,  solarElevation, timestamp); //this triggers Snowpack calculation
Julien Esseiva's avatar
Julien Esseiva committed
182
183
184
185
186
187
188
189
190
		} 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
191
void EnergyBalance::writeSP(const unsigned int max_steps){
192
	if (hasSP()) terrain_radiation->writeSP(max_steps);
193
194
}

195

Julien Esseiva's avatar
Julien Esseiva committed
196
197
198
199
double EnergyBalance::getTiming() const
{
	return timer.getElapsed();
}