WSL/SLF GitLab Repository

EnergyBalance.cc 7.07 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
#include <alpine3d/MPIControl.h>
#include <alpine3d/OMPControl.h>
21

22
23
24
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

32
33
34
35
36
37
38
	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;
39
		OMPControl::getArraySliceParams(nx, nbworkers, ii, thread_startx, thread_nx);
40
		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));
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

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
	}
66
67
}

68

69
70
71
72
EnergyBalance::~EnergyBalance() {
	Destroy( );
}

73
74
75
76
77
78
EnergyBalance& EnergyBalance::operator=(const EnergyBalance& source) {
	if (this != &source) {
		snowpack = source.snowpack;
		terrain_radiation = source.terrain_radiation;
		radfields = source.radfields;
		dem = source.dem;
79
80
		vecMeteo = source.vecMeteo;
		albedo = source.albedo;
81
82
83
		direct = source.direct;
		diffuse = source.diffuse;
		reflected = source.reflected;
84
		direct_unshaded_horizontal = source.reflected;
85
86
87
88
89
90
91
92
		timer = source.timer;
		dimx = source.dimx;
		dimy = source.dimy;
		nbworkers = source.nbworkers;
	}
	return *this;
}

93
94
std::string EnergyBalance::getGridsRequirements() const
{
95
	return "TOP_ALB";
96
97
}

98
99
void EnergyBalance::Destroy()
{
100
101
102
103
	if (terrain_radiation) {
		delete terrain_radiation;
		terrain_radiation = NULL;
	}
104
105
}

Thomas Egger's avatar
Thomas Egger committed
106
void EnergyBalance::setSnowPack(SnowpackInterface& mysnowpack)
107
108
109
110
{
	snowpack = &mysnowpack;
}

Thomas Egger's avatar
Thomas Egger committed
111
112
void EnergyBalance::setAlbedo(const mio::Grid2DObject& in_albedo)
{
113
	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;
119
120
}

121
void EnergyBalance::setStations(const std::vector<mio::MeteoData>& in_vecMeteo)
Thomas Egger's avatar
Thomas Egger committed
122
{
123
	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;
129
130
131
}

void EnergyBalance::setMeteo(const mio::Grid2DObject& in_ilwr,
Thomas Egger's avatar
Thomas Egger committed
132
133
                             const mio::Grid2DObject& in_ta, const mio::Grid2DObject& in_rh, const mio::Grid2DObject& in_p, const mio::Date timestamp)
{
134
	timer.restart();
135

136
	#pragma omp parallel for schedule(dynamic)
137
	for (size_t ii=0; ii<nbworkers; ii++) {
138
		radfields[ii].setStations(vecMeteo, albedo); //calculate the parameters at the radiation stations
139
		size_t startx, nx;
140
141
		radfields[ii].getBandOffsets(startx, nx);
		radfields[ii].setMeteo(mio::Grid2DObject(in_ta, startx, 0, nx, dimy),
142
143
144
		                       mio::Grid2DObject(in_rh, startx, 0, nx, dimy),
		                       mio::Grid2DObject(in_p, startx, 0, nx, dimy),
		                       mio::Grid2DObject(albedo, startx, 0, nx, dimy));
145

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

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

161
	if (terrain_radiation) {
162
		// note: parallelization has to take place inside the TerrainRadiationAlgorithm implementations
163
		terrain_radiation->setMeteo(albedo.grid2D, in_ta.grid2D, in_rh.grid2D, in_ilwr.grid2D);
164
		terrain_radiation->getRadiation(direct, diffuse, reflected, direct_unshaded_horizontal,
165
                                    solarAzimuth, solarElevation);
166
	}
167

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

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

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

177
		timer.stop();
178
		try {
179
			snowpack->setRadiationComponents(global, ilwr, diffuse, reflected,
180
                                       in_ilwr.grid2D,  solarElevation, timestamp); //this triggers Snowpack calculation
181
		} catch(std::exception& e) {
182
			std::cout << "[E] Exception in snowpack->setRadiationComponents()\n";
183
			cout << e.what() << endl;
184
			std::abort(); //force core dump
185
		}
186
187
188
189
	}
	timer.stop();
}

Adrien Michel's avatar
Adrien Michel committed
190
void EnergyBalance::writeSP(const unsigned int max_steps){
191
	if (hasSP()) terrain_radiation->writeSP(max_steps);
192
193
}

194

195
196
197
198
double EnergyBalance::getTiming() const
{
	return timer.getElapsed();
}