Commit edb297e7 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

vert_gen_weights3d1d: OpenMP version.

parent b4c17e7d
Pipeline #4339 failed with stages
in 13 minutes and 16 seconds
......@@ -23,6 +23,7 @@
#include <cdi.h>
#include "cimdOmp.h"
#include "cdo_options.h"
#include "process_int.h"
#include "cdo_vlist.h"
......@@ -104,8 +105,8 @@ vert_interp_lev3d(size_t gridsize, T missval, const Varray<T> &vardata1, Varray<
#endif
for (size_t i = 0; i < gridsize; i++)
{
const auto idx1 = lev_idx1[offset + i];
const auto idx2 = lev_idx2[offset + i];
const auto idx1 = lev_idx1[offset + i] * gridsize + i;
const auto idx2 = lev_idx2[offset + i] * gridsize + i;
const auto wgt1 = lev_wgt1[offset + i];
const auto wgt2 = lev_wgt2[offset + i];
......@@ -165,28 +166,25 @@ vert_gen_weights(int expol, int nlev1, const Varray<double> &lev1, int nlev2, co
if (i1 - 1 == 0) // destination levels is not covert by the first two input z levels
{
lev_idx1[i2] = 1;
lev_idx2[i2] = 1;
lev_idx1[i2] = 0;
lev_idx2[i2] = 0;
lev_wgt1[i2] = 0.0f;
lev_wgt2[i2] = (expol || IS_EQUAL(lev2[i2], val2));
}
else if (i1 == nlev1 - 1) // destination level is beyond the last value of the input z field
{
lev_idx1[i2] = nlev1 - 2;
lev_idx2[i2] = nlev1 - 2;
lev_idx1[i2] = nlev1 - 3;
lev_idx2[i2] = nlev1 - 3;
lev_wgt1[i2] = (expol || IS_EQUAL(lev2[i2], val2));
lev_wgt2[i2] = 0.0f;
}
else // target z values has two bounday values in input z field
{
lev_idx1[i2] = idx1;
lev_idx2[i2] = idx2;
lev_idx1[i2] = idx1 - 1;
lev_idx2[i2] = idx2 - 1;
lev_wgt1[i2] = (lev1[idx2] - lev2[i2]) / (lev1[idx2] - lev1[idx1]);
lev_wgt2[i2] = (lev2[i2] - lev1[idx1]) / (lev1[idx2] - lev1[idx1]);
}
// backshift of the indices because of the two additional levels in input vertical coordinate
lev_idx1[i2]--;
lev_idx2[i2]--;
/*
printf("%d %g %d %d %g %g %d %d %g %g\n", i2, lev2[i2], idx1, idx2, lev1[idx1], lev1[idx2],
lev_idx1[i2], lev_idx2[i2], lev_wgt1[i2], lev_wgt2[i2]);
......@@ -219,30 +217,36 @@ static void
vert_gen_weights3d1d(bool expol, size_t gridsize, int nlev1, const Varray<T> &xlev1, int nlev2, const Varray<double> &lev2,
Varray<int> &xlev_idx1, Varray<int> &xlev_idx2, Varray<float> &xlev_wgt1, Varray<float> &xlev_wgt2)
{
Varray<double> lev1p2(nlev1 + 2);
Varray<float> lev_wgt1(nlev2), lev_wgt2(nlev2);
Varray<int> lev_idx1(nlev2), lev_idx2(nlev2);
const auto nthreads = Threading::ompNumThreads;
Varray2D<double> lev1p2(nthreads, Varray<double>(nlev1 + 2));
Varray2D<float> lev_wgt1(nthreads, Varray<float>(nlev2)), lev_wgt2(nthreads, Varray<float>(nlev2));
Varray2D<int> lev_idx1(nthreads, Varray<int>(nlev2)), lev_idx2(nthreads, Varray<int>(nlev2));
// Check monotony of vertical levels
for (int k = 0; k < nlev1; ++k) lev1p2[k] = xlev1[k * gridsize];
const auto lup = levelDirUp(nlev1, lev1p2.data());
const auto ldown = levelDirDown(nlev1, lev1p2.data());
for (int k = 0; k < nlev1; ++k) lev1p2[0][k] = xlev1[k * gridsize];
const auto lup = levelDirUp(nlev1, lev1p2[0].data());
const auto ldown = levelDirDown(nlev1, lev1p2[0].data());
if (!lup && !ldown) cdoAbort("Non monotonic zaxis!");
const double level_0 = lup ? -1.e33 : 1.e33;
const double level_N = lup ? 1.e33 : -1.e33;
#ifdef _OPENMP
#pragma omp parallel for default(none) shared(expol, gridsize, nlev1, nlev2, lev2, level_0, level_N, lev1p2, xlev1, xlev_idx1, xlev_idx2, xlev_wgt1, xlev_wgt2, lev_idx1, lev_idx2, lev_wgt1, lev_wgt2)
#endif
for (size_t i = 0; i < gridsize; i++)
{
lev1p2[0] = level_0;
lev1p2[nlev1 + 1] = level_N;
for (int k = 0; k < nlev1; ++k) lev1p2[k + 1] = xlev1[k * gridsize + i];
const auto ompthID = cdo_omp_get_thread_num();
lev1p2[ompthID][0] = level_0;
lev1p2[ompthID][nlev1 + 1] = level_N;
for (int k = 0; k < nlev1; ++k) lev1p2[ompthID][k + 1] = xlev1[k * gridsize + i];
vert_gen_weights(expol, nlev1 + 2, lev1p2, nlev2, lev2, lev_idx1, lev_idx2, lev_wgt1, lev_wgt2);
vert_gen_weights(expol, nlev1 + 2, lev1p2[ompthID], nlev2, lev2, lev_idx1[ompthID], lev_idx2[ompthID], lev_wgt1[ompthID], lev_wgt2[ompthID]);
for (int k = 0; k < nlev2; ++k) xlev_idx1[k * gridsize + i] = lev_idx1[k] * gridsize + i;
for (int k = 0; k < nlev2; ++k) xlev_idx2[k * gridsize + i] = lev_idx2[k] * gridsize + i;
for (int k = 0; k < nlev2; ++k) xlev_wgt1[k * gridsize + i] = lev_wgt1[k];
for (int k = 0; k < nlev2; ++k) xlev_wgt2[k * gridsize + i] = lev_wgt2[k];
for (int k = 0; k < nlev2; ++k) xlev_idx1[k * gridsize + i] = lev_idx1[ompthID][k];
for (int k = 0; k < nlev2; ++k) xlev_idx2[k * gridsize + i] = lev_idx2[ompthID][k];
for (int k = 0; k < nlev2; ++k) xlev_wgt1[k * gridsize + i] = lev_wgt1[ompthID][k];
for (int k = 0; k < nlev2; ++k) xlev_wgt2[k * gridsize + i] = lev_wgt2[ompthID][k];
}
}
......
......@@ -70,8 +70,8 @@ vert_gen_weights3d(bool expol, size_t gridsize, int nlev1, Varray<float> &xlev1,
vert_gen_weights(expol, nlev1 + 2, lev1p2, nlev2, lev2, lev_idx1, lev_idx2, lev_wgt1, lev_wgt2);
for (int k = 0; k < nlev2; ++k) xlev_idx1[k * gridsize + i] = lev_idx1[k] * gridsize + i;
for (int k = 0; k < nlev2; ++k) xlev_idx2[k * gridsize + i] = lev_idx2[k] * gridsize + i;
for (int k = 0; k < nlev2; ++k) xlev_idx1[k * gridsize + i] = lev_idx1[k];
for (int k = 0; k < nlev2; ++k) xlev_idx2[k * gridsize + i] = lev_idx2[k];
for (int k = 0; k < nlev2; ++k) xlev_wgt1[k * gridsize + i] = lev_wgt1[k];
for (int k = 0; k < nlev2; ++k) xlev_wgt2[k * gridsize + i] = lev_wgt2[k];
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment