1.3.49
 
Loading...
Searching...
No Matches
EnergyBalanceModel.cu
Go to the documentation of this file.
1
16#include <cuda_runtime.h>
17#include "EnergyBalanceModel.h"
18
19using namespace std;
20
21#define CUDA_CHECK_ERROR(ans) \
22 { \
23 gpuAssert((ans), __FILE__, __LINE__); \
24 }
25inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort = true) {
26 if (code != cudaSuccess) {
27 fprintf(stderr, "GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
28 if (abort)
29 exit(code);
30 }
31}
32
33__device__ float evaluateEnergyBalance(float T, float R, float Qother, float eps, float Ta, float ea, float pressure, float gH, float gS, uint Nsides, float stomatal_sidedness, float heatcapacity, float surfacehumidity, float dt, float Tprev) {
34
35 // Outgoing emission flux
36 float Rout = float(Nsides) * eps * 5.67e-8F * T * T * T * T;
37
38 // Sensible heat flux
39 float QH = cp_air_mol * gH * (T - Ta); // (see Campbell and Norman Eq. 6.8)
40
41 // Latent heat flux
42 float es = 611.0f * expf(17.502f * (T - 273.f) / (T - 273.f + 240.97f));
43 float gM = 1.08f * gH * gS * (stomatal_sidedness / (1.08f * gH + gS * stomatal_sidedness) + (1.f - stomatal_sidedness) / (1.08f * gH + gS * (1.f - stomatal_sidedness)));
44 if (gH == 0 && gS == 0) { // if somehow both go to zero, can get NaN
45 gM = 0;
46 }
47
48 float QL = gM * lambda_mol * (es - ea * surfacehumidity) / pressure;
49
50 // Storage heat flux
51 float storage = 0.f;
52 if (dt > 0) {
53 storage = heatcapacity * (T - Tprev) / dt;
54 }
55
56 // Residual
57 return R - Rout - QH - QL - Qother - storage;
58}
59
60__global__ void solveEnergyBalance(uint Nprimitives, float *To, float *R, float *Qother, float *eps, float *Ta, float *ea, float *pressure, float *gH, float *gS, uint *Nsides, float *stomatal_sidedness, float *TL, float *heatcapacity,
61 float *surfacehumidity, float dt) {
62
63 uint p = blockIdx.x * blockDim.x + threadIdx.x;
64
65 if (p >= Nprimitives) {
66 return;
67 }
68
69 float T;
70
71 float err_max = 0.0001;
72 uint max_iter = 100;
73
74 float T_old_old = To[p];
75
76 float T_old = T_old_old;
77 T_old_old = 400.f;
78
79 float resid_old = evaluateEnergyBalance(T_old, R[p], Qother[p], eps[p], Ta[p], ea[p], pressure[p], gH[p], gS[p], Nsides[p], stomatal_sidedness[p], heatcapacity[p], surfacehumidity[p], dt, To[p]);
80 float resid_old_old = evaluateEnergyBalance(T_old_old, R[p], Qother[p], eps[p], Ta[p], ea[p], pressure[p], gH[p], gS[p], Nsides[p], stomatal_sidedness[p], heatcapacity[p], surfacehumidity[p], dt, To[p]);
81
82 float resid = 100;
83 float err = resid;
84 uint iter = 0;
85 while (err > err_max && iter < max_iter) {
86
87 if (resid_old == resid_old_old) { // this condition will cause NaN
88 err = 0;
89 break;
90 }
91
92 T = fabs((T_old_old * resid_old - T_old * resid_old_old) / (resid_old - resid_old_old));
93
94 resid = evaluateEnergyBalance(T, R[p], Qother[p], eps[p], Ta[p], ea[p], pressure[p], gH[p], gS[p], Nsides[p], stomatal_sidedness[p], heatcapacity[p], surfacehumidity[p], dt, To[p]);
95
96 resid_old_old = resid_old;
97 resid_old = resid;
98
99 // err = fabs(resid);
100 // err = fabs(resid_old-resid_old_old)/fabs(resid_old_old);
101 err = fabs(T_old - T_old_old) / fabs(T_old_old);
102
103 T_old_old = T_old;
104 T_old = T;
105
106 iter++;
107 }
108
109 if (err > err_max) {
110 printf("WARNING (EnergyBalanceModel::solveEnergyBalance): Energy balance did not converge.\n");
111 }
112
113 TL[p] = T;
114}
115
116void EnergyBalanceModel::evaluateSurfaceEnergyBalance(const std::vector<uint> &UUIDs, float dt) {
117
118 //---- Sum up to get total absorbed radiation across all bands ----//
119
120 // Look through all flux primitive data in the context and sum them up in vector Rn. Each element of Rn corresponds to a primitive.
121
122 if (radiation_bands.empty()) {
123 helios::helios_runtime_error("ERROR (EnergyBalanceModel::run): No radiation bands were found.");
124 }
125
126
127 const uint Nprimitives = UUIDs.size();
128
129 std::vector<float> Rn;
130 Rn.resize(Nprimitives, 0);
131
132 std::vector<float> emissivity;
133 emissivity.resize(Nprimitives);
134 for (size_t u = 0; u < Nprimitives; u++) {
135 emissivity.at(u) = 1.f;
136 }
137
138 for (int b = 0; b < radiation_bands.size(); b++) {
139 for (size_t u = 0; u < Nprimitives; u++) {
140 size_t p = UUIDs.at(u);
141
142 char str[50];
143 sprintf(str, "radiation_flux_%s", radiation_bands.at(b).c_str());
144 if (!context->doesPrimitiveDataExist(p, str)) {
145 helios::helios_runtime_error("ERROR (EnergyBalanceModel::run): No radiation was found in the context for band " + std::string(radiation_bands.at(b)) + ". Did you run the radiation model for this band?");
146 } else if (context->getPrimitiveDataType(str) != helios::HELIOS_TYPE_FLOAT) {
147 helios::helios_runtime_error("ERROR (EnergyBalanceModel::run): Radiation primitive data for band " + std::string(radiation_bands.at(b)) + " does not have the correct type of ''float'");
148 }
149 float R;
150 context->getPrimitiveData(p, str, R);
151 Rn.at(u) += R;
152
153 sprintf(str, "emissivity_%s", radiation_bands.at(b).c_str());
154 if (context->doesPrimitiveDataExist(p, str) && context->getPrimitiveDataType(str) == helios::HELIOS_TYPE_FLOAT) {
155 context->getPrimitiveData(p, str, emissivity.at(u));
156 }
157 }
158 }
159
160 //---- Set up temperature solution ----//
161
162 // To,R,Qother,eps,U,L,Ta,ea,pressure,gS,Nsides
163
164 float *To = (float *) malloc(Nprimitives * sizeof(float));
165 float *d_To;
166 CUDA_CHECK_ERROR(cudaMalloc((void **) &d_To, Nprimitives * sizeof(float)));
167
168 float *R = (float *) malloc(Nprimitives * sizeof(float));
169 float *d_R;
170 CUDA_CHECK_ERROR(cudaMalloc((void **) &d_R, Nprimitives * sizeof(float)));
171
172 float *Qother = (float *) malloc(Nprimitives * sizeof(float));
173 float *d_Qother;
174 CUDA_CHECK_ERROR(cudaMalloc((void **) &d_Qother, Nprimitives * sizeof(float)));
175
176 float *eps = (float *) malloc(Nprimitives * sizeof(float));
177 float *d_eps;
178 CUDA_CHECK_ERROR(cudaMalloc((void **) &d_eps, Nprimitives * sizeof(float)));
179
180 float *Ta = (float *) malloc(Nprimitives * sizeof(float));
181 float *d_Ta;
182 CUDA_CHECK_ERROR(cudaMalloc((void **) &d_Ta, Nprimitives * sizeof(float)));
183
184 float *ea = (float *) malloc(Nprimitives * sizeof(float));
185 float *d_ea;
186 CUDA_CHECK_ERROR(cudaMalloc((void **) &d_ea, Nprimitives * sizeof(float)));
187
188 float *pressure = (float *) malloc(Nprimitives * sizeof(float));
189 float *d_pressure;
190 CUDA_CHECK_ERROR(cudaMalloc((void **) &d_pressure, Nprimitives * sizeof(float)));
191
192 float *gH = (float *) malloc(Nprimitives * sizeof(float));
193 float *d_gH;
194 CUDA_CHECK_ERROR(cudaMalloc((void **) &d_gH, Nprimitives * sizeof(float)));
195
196 float *gS = (float *) malloc(Nprimitives * sizeof(float));
197 float *d_gS;
198 CUDA_CHECK_ERROR(cudaMalloc((void **) &d_gS, Nprimitives * sizeof(float)));
199
200 uint *Nsides = (uint *) malloc(Nprimitives * sizeof(uint));
201 uint *d_Nsides;
202 CUDA_CHECK_ERROR(cudaMalloc((void **) &d_Nsides, Nprimitives * sizeof(uint)));
203
204 float *stomatal_sidedness = (float *) malloc(Nprimitives * sizeof(float));
205 float *d_stomatal_sidedness;
206 CUDA_CHECK_ERROR(cudaMalloc((void **) &d_stomatal_sidedness, Nprimitives * sizeof(float)));
207
208 float *heatcapacity = (float *) malloc(Nprimitives * sizeof(float));
209 float *d_heatcapacity;
210 CUDA_CHECK_ERROR(cudaMalloc((void **) &d_heatcapacity, Nprimitives * sizeof(float)));
211
212 float *surfacehumidity = (float *) malloc(Nprimitives * sizeof(float));
213 float *d_surfacehumidity;
214 CUDA_CHECK_ERROR(cudaMalloc((void **) &d_surfacehumidity, Nprimitives * sizeof(float)));
215
216 bool calculated_blconductance_used = false;
217 bool primitive_length_used = false;
218
219 for (uint u = 0; u < Nprimitives; u++) {
220 size_t p = UUIDs.at(u);
221
222 // Initial guess for surface temperature
223 if (context->doesPrimitiveDataExist(p, "temperature") && context->getPrimitiveDataType("temperature") == helios::HELIOS_TYPE_FLOAT) {
224 context->getPrimitiveData(p, "temperature", To[u]);
225 } else {
226 To[u] = temperature_default;
227 }
228 if (To[u] == 0) { // can't have To equal to 0
229 To[u] = 300;
230 }
231
232 // Air temperature
233 if (context->doesPrimitiveDataExist(p, "air_temperature") && context->getPrimitiveDataType("air_temperature") == helios::HELIOS_TYPE_FLOAT) {
234 context->getPrimitiveData(p, "air_temperature", Ta[u]);
235 if (message_flag && Ta[u] < 250.f) {
236 std::cout << "WARNING (EnergyBalanceModel::run): Value of " << Ta[u] << " given in 'air_temperature' primitive data is very small. Values should be given in units of Kelvin. Assuming default value of " << air_temperature_default
237 << std::endl;
238 Ta[u] = air_temperature_default;
239 }
240 } else {
241 Ta[u] = air_temperature_default;
242 }
243
244 // Air relative humidity
245 float hr;
246 if (context->doesPrimitiveDataExist(p, "air_humidity") && context->getPrimitiveDataType("air_humidity") == helios::HELIOS_TYPE_FLOAT) {
247 context->getPrimitiveData(p, "air_humidity", hr);
248 if (hr > 1.f) {
249 if (message_flag) {
250 std::cout << "WARNING (EnergyBalanceModel::run): Value of " << hr << " given in 'air_humidity' primitive data is large than 1. Values should be given as fractional values between 0 and 1. Assuming default value of "
251 << air_humidity_default << std::endl;
252 }
253 hr = air_humidity_default;
254 } else if (hr < 0.f) {
255 if (message_flag) {
256 std::cout << "WARNING (EnergyBalanceModel::run): Value of " << hr << " given in 'air_humidity' primitive data is less than 0. Values should be given as fractional values between 0 and 1. Assuming default value of "
257 << air_humidity_default << std::endl;
258 }
259 hr = air_humidity_default;
260 }
261 } else {
262 hr = air_humidity_default;
263 }
264
265 // Air vapor pressure
266 float esat = esat_Pa(Ta[u]);
267 ea[u] = hr * esat; // Definition of vapor pressure (see Campbell and Norman pp. 42 Eq. 3.11)
268
269 // Air pressure
270 if (context->doesPrimitiveDataExist(p, "air_pressure") && context->getPrimitiveDataType("air_pressure") == helios::HELIOS_TYPE_FLOAT) {
271 context->getPrimitiveData(p, "air_pressure", pressure[u]);
272 if (pressure[u] < 10000.f) {
273 if (message_flag) {
274 std::cout << "WARNING (EnergyBalanceModel::run): Value of " << pressure[u] << " given in 'air_pressure' primitive data is very small. Values should be given in units of Pascals. Assuming default value of " << pressure_default
275 << std::endl;
276 }
277 pressure[u] = pressure_default;
278 }
279 } else {
280 pressure[u] = pressure_default;
281 }
282
283 // Number of sides emitting radiation
284 Nsides[u] = 2; // default is 2
285 if (context->doesPrimitiveDataExist(p, "twosided_flag") && context->getPrimitiveDataType("twosided_flag") == helios::HELIOS_TYPE_UINT) {
286 uint flag;
287 context->getPrimitiveData(p, "twosided_flag", flag);
288 if (flag == 0) {
289 Nsides[u] = 1;
290 }
291 }
292
293 // Number of evaporating/transpiring faces
294 stomatal_sidedness[u] = 0.f; // if Nsides=1, force this to be 0 (all stomata on upper surface)
295 if (Nsides[u] == 2 && context->doesPrimitiveDataExist(p, "stomatal_sidedness") && context->getPrimitiveDataType("stomatal_sidedness") == helios::HELIOS_TYPE_FLOAT) {
296 context->getPrimitiveData(p, "stomatal_sidedness", stomatal_sidedness[u]);
297 // this is for backward compatability prior to v1.3.17
298 } else if (Nsides[u] == 2 && context->doesPrimitiveDataExist(p, "evaporating_faces") && context->getPrimitiveDataType("evaporating_faces") == helios::HELIOS_TYPE_UINT) {
299 uint flag;
300 context->getPrimitiveData(p, "evaporating_faces", flag);
301 if (flag == 1) { // stomata on one side
302 stomatal_sidedness[u] = 0.f;
303 } else if (flag == 2) {
304 stomatal_sidedness[u] = 0.5f;
305 }
306 }
307
308 // Boundary-layer conductance to heat
309 if (context->doesPrimitiveDataExist(p, "boundarylayer_conductance") && context->getPrimitiveDataType("boundarylayer_conductance") == helios::HELIOS_TYPE_FLOAT) {
310 context->getPrimitiveData(p, "boundarylayer_conductance", gH[u]);
311 } else {
312
313 // Wind speed
314 float U;
315 if (context->doesPrimitiveDataExist(p, "wind_speed") && context->getPrimitiveDataType("wind_speed") == helios::HELIOS_TYPE_FLOAT) {
316 context->getPrimitiveData(p, "wind_speed", U);
317 } else {
318 U = wind_speed_default;
319 }
320
321 // Characteristic size of primitive
322 float L;
323 if (context->doesPrimitiveDataExist(p, "object_length") && context->getPrimitiveDataType("object_length") == helios::HELIOS_TYPE_FLOAT) {
324 context->getPrimitiveData(p, "object_length", L);
325 if (L == 0) {
326 L = sqrt(context->getPrimitiveArea(p));
327 primitive_length_used = true;
328 }
329 } else if (context->getPrimitiveParentObjectID(p) > 0) {
330 uint objID = context->getPrimitiveParentObjectID(p);
331 L = sqrt(context->getObjectArea(objID));
332 } else {
333 L = sqrt(context->getPrimitiveArea(p));
334 primitive_length_used = true;
335 }
336
337 gH[u] = 0.135f * sqrt(U / L) * float(Nsides[u]);
338
339 calculated_blconductance_used = true;
340 }
341
342 // Moisture conductance
343 if (context->doesPrimitiveDataExist(p, "moisture_conductance") && context->getPrimitiveDataType("moisture_conductance") == helios::HELIOS_TYPE_FLOAT) {
344 context->getPrimitiveData(p, "moisture_conductance", gS[u]);
345 } else {
346 gS[u] = gS_default;
347 }
348
349 // Other fluxes
350 if (context->doesPrimitiveDataExist(p, "other_surface_flux") && context->getPrimitiveDataType("other_surface_flux") == helios::HELIOS_TYPE_FLOAT) {
351 context->getPrimitiveData(p, "other_surface_flux", Qother[u]);
352 } else {
353 Qother[u] = Qother_default;
354 }
355
356 // Object heat capacity
357 if (context->doesPrimitiveDataExist(p, "heat_capacity") && context->getPrimitiveDataType("heat_capacity") == helios::HELIOS_TYPE_FLOAT) {
358 context->getPrimitiveData(p, "heat_capacity", heatcapacity[u]);
359 } else {
360 heatcapacity[u] = heatcapacity_default;
361 }
362
363 // Surface humidity
364 if (context->doesPrimitiveDataExist(p, "surface_humidity") && context->getPrimitiveDataType("surface_humidity") == helios::HELIOS_TYPE_FLOAT) {
365 context->getPrimitiveData(p, "surface_humidity", surfacehumidity[u]);
366 } else {
367 surfacehumidity[u] = surface_humidity_default;
368 }
369
370 // Emissivity
371 eps[u] = emissivity.at(u);
372
373 // Net absorbed radiation
374 R[u] = Rn.at(u);
375 }
376
377 // if we used the calculated boundary-layer conductance, enable output primitive data "boundarylayer_conductance_out" so that it can be used by other plug-ins
378 if (calculated_blconductance_used) {
379 auto it = find(output_prim_data.begin(), output_prim_data.end(), "boundarylayer_conductance_out");
380 if (it == output_prim_data.end()) {
381 output_prim_data.emplace_back("boundarylayer_conductance_out");
382 }
383 }
384
385 // if the length of a primitive that is not a member of an object was used, issue a warning
386 if (message_flag && primitive_length_used) {
387 std::cout
388 << "WARNING (EnergyBalanceModel::run): The length of a primitive that is not a member of a compound object was used to calculate the boundary-layer conductance. This often results in incorrect values because the length should be that of the object (e.g., leaf, stem) not the primitive. Make sure this is what you intended."
389 << std::endl;
390 }
391
392 // To,R,Qother,eps,U,L,Ta,ea,pressure,gS,Nsides
393 CUDA_CHECK_ERROR(cudaMemcpy(d_To, To, Nprimitives * sizeof(float), cudaMemcpyHostToDevice));
394 CUDA_CHECK_ERROR(cudaMemcpy(d_R, R, Nprimitives * sizeof(float), cudaMemcpyHostToDevice));
395 CUDA_CHECK_ERROR(cudaMemcpy(d_Qother, Qother, Nprimitives * sizeof(float), cudaMemcpyHostToDevice));
396 CUDA_CHECK_ERROR(cudaMemcpy(d_eps, eps, Nprimitives * sizeof(float), cudaMemcpyHostToDevice));
397 CUDA_CHECK_ERROR(cudaMemcpy(d_Ta, Ta, Nprimitives * sizeof(float), cudaMemcpyHostToDevice));
398 CUDA_CHECK_ERROR(cudaMemcpy(d_ea, ea, Nprimitives * sizeof(float), cudaMemcpyHostToDevice));
399 CUDA_CHECK_ERROR(cudaMemcpy(d_pressure, pressure, Nprimitives * sizeof(float), cudaMemcpyHostToDevice));
400 CUDA_CHECK_ERROR(cudaMemcpy(d_gH, gH, Nprimitives * sizeof(float), cudaMemcpyHostToDevice));
401 CUDA_CHECK_ERROR(cudaMemcpy(d_gS, gS, Nprimitives * sizeof(float), cudaMemcpyHostToDevice));
402 CUDA_CHECK_ERROR(cudaMemcpy(d_Nsides, Nsides, Nprimitives * sizeof(uint), cudaMemcpyHostToDevice));
403 CUDA_CHECK_ERROR(cudaMemcpy(d_stomatal_sidedness, stomatal_sidedness, Nprimitives * sizeof(float), cudaMemcpyHostToDevice));
404 CUDA_CHECK_ERROR(cudaMemcpy(d_heatcapacity, heatcapacity, Nprimitives * sizeof(float), cudaMemcpyHostToDevice));
405 CUDA_CHECK_ERROR(cudaMemcpy(d_surfacehumidity, surfacehumidity, Nprimitives * sizeof(float), cudaMemcpyHostToDevice));
406
407 float *T = (float *) malloc(Nprimitives * sizeof(float));
408 float *d_T;
409 CUDA_CHECK_ERROR(cudaMalloc((void **) &d_T, Nprimitives * sizeof(float)));
410
411 // launch kernel
412 dim3 dimBlock(64, 1);
413 dim3 dimGrid(ceil(Nprimitives / 64.f));
414
415 solveEnergyBalance<<<dimGrid, dimBlock>>>(Nprimitives, d_To, d_R, d_Qother, d_eps, d_Ta, d_ea, d_pressure, d_gH, d_gS, d_Nsides, d_stomatal_sidedness, d_T, d_heatcapacity, d_surfacehumidity, dt);
416
417 CUDA_CHECK_ERROR(cudaPeekAtLastError());
418 CUDA_CHECK_ERROR(cudaDeviceSynchronize());
419
420 CUDA_CHECK_ERROR(cudaMemcpy(T, d_T, Nprimitives * sizeof(float), cudaMemcpyDeviceToHost));
421
422 for (uint u = 0; u < Nprimitives; u++) {
423 size_t UUID = UUIDs.at(u);
424
425 if (T[u] != T[u]) {
426 T[u] = temperature_default;
427 }
428
429 context->setPrimitiveData(UUID, "temperature", T[u]);
430
431 float QH = cp_air_mol * gH[u] * (T[u] - Ta[u]);
432 context->setPrimitiveData(UUID, "sensible_flux", QH);
433
434 float es = esat_Pa(T[u]);
435 float gM = 1.08f * gH[u] * gS[u] * (stomatal_sidedness[u] / (1.08f * gH[u] + gS[u] * stomatal_sidedness[u]) + (1.f - stomatal_sidedness[u]) / (1.08f * gH[u] + gS[u] * (1.f - stomatal_sidedness[u])));
436 if (gH[u] == 0 && gS[u] == 0) { // if somehow both go to zero, can get NaN
437 gM = 0;
438 }
439 float QL = lambda_mol * gM * (es - ea[u]) / pressure[u];
440 context->setPrimitiveData(UUID, "latent_flux", QL);
441
442 for (const auto &data_label: output_prim_data) {
443 if (data_label == "boundarylayer_conductance_out") {
444 context->setPrimitiveData(UUID, "boundarylayer_conductance_out", gH[u]);
445 } else if (data_label == "vapor_pressure_deficit") {
446 float vpd = (es - ea[u]) / pressure[u];
447 context->setPrimitiveData(UUID, "vapor_pressure_deficit", vpd);
448 } else if (data_label == "net_radiation_flux") {
449 float Rnet = R[u] - float(Nsides[u]) * eps[u] * 5.67e-8F * std::pow(T[u], 4);
450 context->setPrimitiveData(UUID, "net_radiation_flux", Rnet);
451 } else if (data_label == "storage_flux") {
452 float storage = 0.f;
453 if (dt > 0) {
454 storage = heatcapacity[u] * (T[u] - To[u]) / dt;
455 }
456 context->setPrimitiveData(UUID, "storage_flux", storage);
457 }
458 }
459 }
460
461 free(To);
462 free(R);
463 free(Qother);
464 free(eps);
465 free(Ta);
466 free(ea);
467 free(pressure);
468 free(gH);
469 free(gS);
470 free(Nsides);
471 free(stomatal_sidedness);
472 free(heatcapacity);
473 free(surfacehumidity);
474 free(T);
475
476
477 CUDA_CHECK_ERROR(cudaFree(d_To));
478 CUDA_CHECK_ERROR(cudaFree(d_R));
479 CUDA_CHECK_ERROR(cudaFree(d_Qother));
480 CUDA_CHECK_ERROR(cudaFree(d_eps));
481 CUDA_CHECK_ERROR(cudaFree(d_Ta));
482 CUDA_CHECK_ERROR(cudaFree(d_ea));
483 CUDA_CHECK_ERROR(cudaFree(d_pressure));
484 CUDA_CHECK_ERROR(cudaFree(d_gH));
485 CUDA_CHECK_ERROR(cudaFree(d_gS));
486 CUDA_CHECK_ERROR(cudaFree(d_Nsides));
487 CUDA_CHECK_ERROR(cudaFree(d_stomatal_sidedness));
488 CUDA_CHECK_ERROR(cudaFree(d_heatcapacity));
489 CUDA_CHECK_ERROR(cudaFree(d_surfacehumidity));
490 CUDA_CHECK_ERROR(cudaFree(d_T));
491}