diff --git a/CMakeLists.txt b/CMakeLists.txt index 37e4b6016b..d50966f5e4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -185,7 +185,7 @@ adios_option(Sodium "Enable support for Sodium for encryption" AUTO) adios_option(Catalyst "Enable support for in situ visualization plugin using ParaView Catalyst" AUTO) adios_option(Campaign "Enable support for Campaigns (requires SQLite3 and ZLIB)" AUTO) adios_option(AWSSDK "Enable support for S3 compatible storage using AWS SDK's S3 module" OFF) -adios_option(Derived_Variable "Enable support for derived variables" OFF) +adios_option(Derived_Variable "Enable support for derived variables" ON) adios_option(PIP "Enable support for pip packaging" OFF) option(ADIOS2_Blosc2_PREFER_SHARED "Prefer shared Blosc2 libraries" ON) diff --git a/source/adios2/toolkit/derived/ExprHelper.h b/source/adios2/toolkit/derived/ExprHelper.h index e1d31b8bff..38672cfec4 100644 --- a/source/adios2/toolkit/derived/ExprHelper.h +++ b/source/adios2/toolkit/derived/ExprHelper.h @@ -21,7 +21,8 @@ enum ExpressionOperator OP_ADD, OP_SQRT, OP_POW, - OP_MAGN + OP_MAGN, + OP_CURL }; struct OperatorProperty @@ -39,6 +40,7 @@ const std::map op_property = { {ExpressionOperator::OP_ADD, {"ADD", true}}, {ExpressionOperator::OP_SQRT, {"SQRT", false}}, {ExpressionOperator::OP_POW, {"POW", false}}, + {ExpressionOperator::OP_CURL, {"CURL", false}}, {ExpressionOperator::OP_MAGN, {"MAGNITUDE", false}}}; const std::map string_to_op = { @@ -49,6 +51,7 @@ const std::map string_to_op = { {"add", ExpressionOperator::OP_ADD}, {"ADD", ExpressionOperator::OP_ADD}, {"SQRT", ExpressionOperator::OP_SQRT}, {"sqrt", ExpressionOperator::OP_SQRT}, {"POW", ExpressionOperator::OP_POW}, {"^", ExpressionOperator::OP_POW}, + {"CURL", ExpressionOperator::OP_CURL}, {"curl", ExpressionOperator::OP_CURL}, {"MAGNITUDE", ExpressionOperator::OP_MAGN}, {"magnitude", ExpressionOperator::OP_MAGN}}; inline std::string get_op_name(ExpressionOperator op) { return op_property.at(op).name; } diff --git a/source/adios2/toolkit/derived/Function.cpp b/source/adios2/toolkit/derived/Function.cpp index cb145d8927..441be6a7c6 100644 --- a/source/adios2/toolkit/derived/Function.cpp +++ b/source/adios2/toolkit/derived/Function.cpp @@ -54,15 +54,119 @@ Dims SameDimsFunc(std::vector input) // check that all dimenstions are the same if (input.size() > 1) { - bool dim_are_equal = std::equal(input.begin() + 1, input.end(), input.begin()); + Dims first_element = input[0]; + bool dim_are_equal = std::all_of(input.begin() + 1, input.end(), + [&first_element](Dims x) { return x == first_element; }); if (!dim_are_equal) - helper::Throw("Derived", "Function", "SameDimFunc", + helper::Throw("Derived", "Function", "SameDimsFunc", "Invalid variable dimensions"); } // return the first dimension return input[0]; } +// Input Dims are the same, output is combination of all inputs +Dims CurlDimsFunc(std::vector input) +{ + // check that all dimenstions are the same + if (input.size() > 1) + { + Dims first_element = input[0]; + bool dim_are_equal = std::all_of(input.begin() + 1, input.end(), + [&first_element](Dims x) { return x == first_element; }); + if (!dim_are_equal) + helper::Throw("Derived", "Function", "CurlDimsFunc", + "Invalid variable dimensions"); + } + // return original dimensions with added dimension of number of inputs + Dims output = input[0]; + output.push_back(input.size()); + return output; +} + +/* + * Linear Interpolation - average difference around point "index" + * can be used to approximate derivatives + * + * Input: + * input - data assumed to be uniform/densely populated + * index - index of point of interest + * dim - in which dimension we are approximating the partial derivative + */ +// template +float linear_interp(DerivedData input, size_t index, size_t dim) +{ + size_t stride = 1; + size_t range; + size_t offset; + float result; + float *data = (float *)input.Data; + + for (size_t i = 0; i < input.Count.size() - (dim + 1); ++i) + { + stride *= input.Count[input.Count.size() - (i + 1)]; + } + size_t ind1 = index - stride; + size_t ind2 = index + stride; + range = stride * input.Count[dim]; + offset = index % range; + + if ((offset < stride) && (range - offset <= stride)) + { + return 0; + } + else if (offset < stride) + { + result = data[ind2] - data[index]; + } + else if (range - offset <= stride) + { + result = data[index] - data[ind1]; + } + else + { + result = (data[ind2] - data[ind1]) / 2; + } + + return result; +} + +/* + * Input: 3D vector field F(x,y,z)= {F1(x,y,z), F2(x,y,z), F3(x,y,z)} + * + * inputData - (3) components of 3D vector field + * + * Computation: + * curl(F(x,y,z)) = (partial(F3,y) - partial(F2,z))i + * + (partial(F1,z) - partial(F3,x))j + * + (partial(F2,x) - partial(F1,y))k + * + * boundaries are calculated only with data in block + * (ex: partial derivatives in x direction at point (0,0,0) + * only use data from (1,0,0), etc ) + */ +DerivedData Curl3DFunc(const std::vector inputData, DataType type) +{ + size_t dataSize = inputData[0].Count[0] * inputData[0].Count[1] * inputData[0].Count[2]; + + DerivedData curl; + // ToDo - template type + float *data = (float *)malloc(dataSize * sizeof(float) * 3); + curl.Start = inputData[0].Start; + curl.Start.push_back(0); + curl.Count = inputData[0].Count; + curl.Count.push_back(3); + + for (size_t i = 0; i < dataSize; ++i) + { + data[3 * i] = linear_interp(inputData[2], i, 1) - linear_interp(inputData[1], i, 2); + data[3 * i + 1] = linear_interp(inputData[0], i, 2) - linear_interp(inputData[2], i, 0); + data[3 * i + 2] = linear_interp(inputData[1], i, 0) - linear_interp(inputData[0], i, 1); + } + curl.Data = data; + return curl; +} + #define declare_template_instantiation(T) \ T *ApplyOneToOne(std::vector, size_t, std::function); diff --git a/source/adios2/toolkit/derived/Function.h b/source/adios2/toolkit/derived/Function.h index 5dfa5aba97..d41aab6581 100644 --- a/source/adios2/toolkit/derived/Function.h +++ b/source/adios2/toolkit/derived/Function.h @@ -26,11 +26,17 @@ struct OperatorFunctions DerivedData AddFunc(std::vector input, DataType type); DerivedData MagnitudeFunc(std::vector input, DataType type); +DerivedData Curl3DFunc(std::vector input, DataType type); + +template +T linear_interp(T *data, size_t index, size_t count, size_t stride = 1); Dims SameDimsFunc(std::vector input); +Dims CurlDimsFunc(std::vector input); const std::map OpFunctions = { {adios2::detail::ExpressionOperator::OP_ADD, {AddFunc, SameDimsFunc}}, + {adios2::detail::ExpressionOperator::OP_CURL, {Curl3DFunc, CurlDimsFunc}}, {adios2::detail::ExpressionOperator::OP_MAGN, {MagnitudeFunc, SameDimsFunc}}}; template diff --git a/source/adios2/toolkit/derived/parser/pregen-source/lexer.cpp b/source/adios2/toolkit/derived/parser/pregen-source/lexer.cpp index 11c2e6b66a..598731e63c 100644 --- a/source/adios2/toolkit/derived/parser/pregen-source/lexer.cpp +++ b/source/adios2/toolkit/derived/parser/pregen-source/lexer.cpp @@ -1903,4 +1903,3 @@ adios2::detail::ASTDriver::parse (const std::string input) parse.set_debug_level (trace_parsing); parse (); } - diff --git a/testing/adios2/derived/TestBPDerivedCorrectness.cpp b/testing/adios2/derived/TestBPDerivedCorrectness.cpp index 886b942db3..9d96cc5f41 100644 --- a/testing/adios2/derived/TestBPDerivedCorrectness.cpp +++ b/testing/adios2/derived/TestBPDerivedCorrectness.cpp @@ -3,7 +3,6 @@ #include #include -#include #include #include #include @@ -169,6 +168,157 @@ TEST(DerivedCorrectness, MagCorrectnessTest) bpFileReader.Close(); } +TEST(DerivedCorrectness, CurlCorrectnessTest) +{ + const size_t Nx = 25, Ny = 70, Nz = 13; + float error_limit = 0.0000001f; + + // Application variable + std::vector simArray1(Nx * Ny * Nz); + std::vector simArray2(Nx * Ny * Nz); + std::vector simArray3(Nx * Ny * Nz); + for (size_t i = 0; i < Nx; ++i) + { + for (size_t j = 0; j < Ny; ++j) + { + for (size_t k = 0; k < Nz; ++k) + { + size_t idx = (i * Ny * Nz) + (j * Nz) + k; + float x = static_cast(i); + float y = static_cast(j); + float z = static_cast(k); + // Linear curl example + simArray1[idx] = (6 * x * y) + (7 * z); + simArray2[idx] = (4 * x * z) + powf(y, 2); + simArray3[idx] = sqrtf(z) + (2 * x * y); + /* Less linear example + simArray1[idx] = sinf(z); + simArray2[idx] = 4 * x; + simArray3[idx] = powf(y, 2) * cosf(x); + */ + /* Nonlinear example + simArray1[idx] = expf(2 * y) * sinf(x); + simArray2[idx] = sqrtf(z + 1) * cosf(x); + simArray3[idx] = powf(x, 2) * sinf(y) + (6 * z); + */ + } + } + } + + adios2::ADIOS adios; + adios2::IO bpOut = adios.DeclareIO("BPWriteExpression"); + std::vector varname = {"sim3/VX", "sim3/VY", "sim3/VZ"}; + std::string derivedname = "derived/curlV"; + + auto VX = bpOut.DefineVariable(varname[0], {Nx, Ny, Nz}, {0, 0, 0}, {Nx, Ny, Nz}); + auto VY = bpOut.DefineVariable(varname[1], {Nx, Ny, Nz}, {0, 0, 0}, {Nx, Ny, Nz}); + auto VZ = bpOut.DefineVariable(varname[2], {Nx, Ny, Nz}, {0, 0, 0}, {Nx, Ny, Nz}); + // clang-format off + bpOut.DefineDerivedVariable(derivedname, + "Vx =" + varname[0] + " \n" + "Vy =" + varname[1] + " \n" + "Vz =" + varname[2] + " \n" + "curl(Vx,Vy,Vz)", + adios2::DerivedVarType::StoreData); + // clang-format on + std::string filename = "expCurl.bp"; + adios2::Engine bpFileWriter = bpOut.Open(filename, adios2::Mode::Write); + + bpFileWriter.BeginStep(); + bpFileWriter.Put(VX, simArray1.data()); + bpFileWriter.Put(VY, simArray2.data()); + bpFileWriter.Put(VZ, simArray3.data()); + bpFileWriter.EndStep(); + bpFileWriter.Close(); + + adios2::IO bpIn = adios.DeclareIO("BPReadCurlExpression"); + adios2::Engine bpFileReader = bpIn.Open(filename, adios2::Mode::Read); + + std::vector readVX; + std::vector readVY; + std::vector readVZ; + // TODO/DEBUG - VERIFY DATATYPE + std::vector readCurl; + + std::vector> calcCurl; + double sum_x = 0; + double sum_y = 0; + double sum_z = 0; + bpFileReader.BeginStep(); + auto varVX = bpIn.InquireVariable(varname[0]); + auto varVY = bpIn.InquireVariable(varname[1]); + auto varVZ = bpIn.InquireVariable(varname[2]); + auto varCurl = bpIn.InquireVariable(derivedname); + + bpFileReader.Get(varVX, readVX); + bpFileReader.Get(varVY, readVY); + bpFileReader.Get(varVZ, readVZ); + bpFileReader.Get(varCurl, readCurl); + bpFileReader.EndStep(); + + float curl_x, curl_y, curl_z; + float err_x, err_y, err_z; + for (size_t i = 0; i < Nx; ++i) + { + for (size_t j = 0; j < Ny; ++j) + { + for (size_t k = 0; k < Nz; ++k) + { + size_t idx = (i * Ny * Nz) + (j * Nz) + k; + float x = static_cast(i); + float y = static_cast(j); + float z = static_cast(k); + // Linear example + curl_x = -(2 * x); + curl_y = 7 - (2 * y); + curl_z = (4 * z) - (6 * x); + /* Less linear + curl_x = 2 * y * cosf(x); + curl_y = cosf(z) + (powf(y, 2) * sinf(x)); + curl_z = 4; + */ + /* Nonlinear example + curl_x = powf(x, 2) * cosf(y) - (cosf(x) / (2 * sqrtf(z + 1))); + curl_y = -2 * x * sinf(y); + curl_z = -sqrtf(z + 1) * sinf(x) - (2 * expf(2 * y) * sinf(x)); + */ + if (fabs(curl_x) < 1) + { + err_x = fabs(curl_x - readCurl[3 * idx]) / (1 + fabs(curl_x)); + } + else + { + err_x = fabs(curl_x - readCurl[3 * idx]) / fabs(curl_x); + } + if (fabs(curl_y) < 1) + { + err_y = fabs(curl_y - readCurl[3 * idx + 1]) / (1 + fabs(curl_y)); + } + else + { + err_y = fabs(curl_y - readCurl[3 * idx + 1]) / fabs(curl_y); + } + if (fabs(curl_z) < 1) + { + err_z = fabs(curl_z - readCurl[3 * idx + 2]) / (1 + fabs(curl_z)); + } + else + { + err_z = fabs(curl_z - readCurl[3 * idx + 2]) / fabs(curl_z); + } + sum_x += err_x; + sum_y += err_y; + sum_z += err_z; + } + } + } + bpFileReader.Close(); + EXPECT_LT((sum_x + sum_y + sum_z) / (3 * Nx * Ny * Nz), error_limit); + EXPECT_LT(sum_x / (Nx * Ny * Nz), error_limit); + EXPECT_LT(sum_y / (Nx * Ny * Nz), error_limit); + EXPECT_LT(sum_z / (Nx * Ny * Nz), error_limit); +} + int main(int argc, char **argv) { int result;