WSL/SLF GitLab Repository

EnergyBalance.cc 7.32 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
29
30
              : snowpack(NULL), terrain_radiation(NULL), radfields(), dem(dem_in), vecMeteo(), 
                albedo(dem, 0.), direct_unshaded_horizontal(dem_in.getNx(), dem_in.getNy(), 0.),
                direct(dem_in.getNx(), dem_in.getNy(), 0.), diffuse(dem_in.getNx(), dem_in.getNy(), 0.), 
                reflected(dem_in.getNx(), dem_in.getNy(), 0.), timer(), cfg(cfg_in), 
                dimx(dem_in.getNx()), dimy(dem_in.getNy()), nbworkers(i_nbworkers)
31
{
32

Julien Esseiva's avatar
Julien Esseiva committed
33
34
35
36
37
38
39
40
41
	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;
42
43
		std::cout << "[i] EnergyBalance worker " << ii << " on process " << instance.rank() << " will start at offset " <<
		offset << " with nx " << thread_nx << "\n";
44
		radfields.push_back(RadiationField(dem_in, offset, thread_nx));
Julien Esseiva's avatar
Julien Esseiva committed
45
46
	}

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

Julien Esseiva's avatar
Julien Esseiva committed
50
	// Every MPI process will have its own copy of terrain_radiation object with full DEM
51
52
53
54
55
	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";

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

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

void EnergyBalance::setMeteo(const mio::Grid2DObject& in_ilwr,
133
134
                             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
135
136
137
138
139
{
	timer.restart();

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

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

160
161
162
163
	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
164
165
	if (terrain_radiation) {
		// note: parallelization has to take place inside the TerrainRadiationAlgorithm implementations
166
		terrain_radiation->setMeteo(albedo.grid2D, in_ta.grid2D);
167
		terrain_radiation->getRadiation(direct, diffuse, reflected, direct_unshaded_horizontal,
168
                                    in_ilwr.grid2D,sky_ilwr,terrain_ilwr, solarAzimuth, solarElevation);
169
170
171
		if (hasSP()){
			terrain_radiation->setSP(radfields[0].getDate(), solarAzimuth, solarElevation);
		}
Julien Esseiva's avatar
Julien Esseiva committed
172
173
174
175
176
177
178
179
180
181
182
183
184
	}

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

200

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