WSL/SLF GitLab Repository

EnergyBalance.cc 7.23 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
	// Every MPI process will have its own copy of terrain_radiation object with full DEM
50
51
52
53
54
	terrain_radiation = TerrainRadiationFactory::getAlgorithm(cfg, dem, nbworkers);
	const std::string algo = terrain_radiation->algo;
	if (instance.master())
		std::cout << "[i] Using terrain radiation with model: " << algo << "\n";

55
56
57
	bool write_sky_vf=false;
	cfg.getValue("WRITE_SKY_VIEW_FACTOR", "output", write_sky_vf,IOUtils::nothrow);

Adrien Michel's avatar
Adrien Michel committed
58
	if(write_sky_vf){
59
60
		std::cout << "[i] Writing sky view factor grid" << std::endl;
		mio::IOManager io(cfg);
Adrien Michel's avatar
Adrien Michel committed
61
		mio::Array2D<double> sky_vf(dimx,dimy,0);
62
		terrain_radiation->getSkyViewFactor(sky_vf);
Adrien Michel's avatar
Adrien Michel committed
63
64
		if(MPIControl::instance().master())
			io.write2DGrid(mio::Grid2DObject(dem_in.cellsize,dem_in.llcorner,sky_vf), "SKY_VIEW_FACTOR");
65
	}
Julien Esseiva's avatar
Julien Esseiva committed
66
67
}

68

Julien Esseiva's avatar
Julien Esseiva committed
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
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;
84
		direct_unshaded_horizontal = source.reflected;
Julien Esseiva's avatar
Julien Esseiva committed
85
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
		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;
114
115
116
117
118
	//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
119
120
121
122
123
}

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

void EnergyBalance::setMeteo(const mio::Grid2DObject& in_ilwr,
132
133
                             const mio::Grid2DObject& in_ta, const mio::Grid2DObject& in_rh,
                             const mio::Grid2DObject& in_p, const mio::Date timestamp)
Julien Esseiva's avatar
Julien Esseiva committed
134
135
136
137
138
{
	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);

162
163
164
165
	mio::Array2D<double> sky_ilwr(in_ilwr.grid2D);
	mio::Array2D<double> terrain_ilwr(in_ilwr.grid2D);
	sky_ilwr=0;
	terrain_ilwr=0;
Julien Esseiva's avatar
Julien Esseiva committed
166
167
	if (terrain_radiation) {
		// note: parallelization has to take place inside the TerrainRadiationAlgorithm implementations
168
		terrain_radiation->setMeteo(albedo.grid2D, in_ta.grid2D);
169
		terrain_radiation->getRadiation(direct, diffuse, reflected, direct_unshaded_horizontal,
170
                                    in_ilwr.grid2D,sky_ilwr,terrain_ilwr, solarAzimuth, solarElevation);
Julien Esseiva's avatar
Julien Esseiva committed
171
172
173
174
175
176
177
178
179
180
181
182
183
	}

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

199

Julien Esseiva's avatar
Julien Esseiva committed
200
201
202
203
double EnergyBalance::getTiming() const
{
	return timer.getElapsed();
}