Add the CURL function to derived variables

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)

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<ExpressionOperator, OperatorProperty> 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<std::string, ExpressionOperator> string_to_op = {
@@ -49,6 +51,7 @@ const std::map<std::string, ExpressionOperator> 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; }

source/adios2/toolkit/derived/Function.cpp

@@ -63,6 +63,106 @@ Dims SameDimsFunc(std::vector<Dims> input)
     return input[0];
 }
 
+// Input Dims are the same, output is combination of all inputs
+Dims AggrSameDimsFunc(std::vector<Dims> 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());
+        if (!dim_are_equal)
+            helper::Throw<std::invalid_argument>("Derived", "Function", "AggrSameDimFunc",
+                                                 "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 <class T>
+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<DerivedData> 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<DerivedData>, size_t, std::function<T(T, T)>);
 

source/adios2/toolkit/derived/Function.h

@@ -26,11 +26,17 @@ struct OperatorFunctions
 
 DerivedData AddFunc(std::vector<DerivedData> input, DataType type);
 DerivedData MagnitudeFunc(std::vector<DerivedData> input, DataType type);
+DerivedData Curl3DFunc(std::vector<DerivedData> input, DataType type);
+
+template <class T>
+T linear_interp(T *data, size_t index, size_t count, size_t stride = 1);
 
 Dims SameDimsFunc(std::vector<Dims> input);
+Dims AggrSameDimsFunc(std::vector<Dims> input);
 
 const std::map<adios2::detail::ExpressionOperator, OperatorFunctions> OpFunctions = {
     {adios2::detail::ExpressionOperator::OP_ADD, {AddFunc, SameDimsFunc}},
+    {adios2::detail::ExpressionOperator::OP_CURL, {Curl3DFunc, AggrSameDimsFunc}},
     {adios2::detail::ExpressionOperator::OP_MAGN, {MagnitudeFunc, SameDimsFunc}}};
 
 template <class T>

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 ();
 }
-

testing/adios2/derived/TestBPDerivedCorrectness.cpp

@@ -3,7 +3,6 @@
 
 #include <cmath>
 #include <iostream>
-#include <limits>
 #include <numeric>
 #include <random>
 #include <stdexcept>
@@ -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<float> simArray1(Nx * Ny * Nz);
+    std::vector<float> simArray2(Nx * Ny * Nz);
+    std::vector<float> 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<float>(i);
+                float y = static_cast<float>(j);
+                float z = static_cast<float>(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<std::string> varname = {"sim3/VX", "sim3/VY", "sim3/VZ"};
+    std::string derivedname = "derived/curlV";
+
+    auto VX = bpOut.DefineVariable<float>(varname[0], {Nx, Ny, Nz}, {0, 0, 0}, {Nx, Ny, Nz});
+    auto VY = bpOut.DefineVariable<float>(varname[1], {Nx, Ny, Nz}, {0, 0, 0}, {Nx, Ny, Nz});
+    auto VZ = bpOut.DefineVariable<float>(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<float> readVX;
+    std::vector<float> readVY;
+    std::vector<float> readVZ;
+    // TODO/DEBUG - VERIFY DATATYPE
+    std::vector<float> readCurl;
+
+    std::vector<std::vector<float>> calcCurl;
+    double sum_x = 0;
+    double sum_y = 0;
+    double sum_z = 0;
+    bpFileReader.BeginStep();
+    auto varVX = bpIn.InquireVariable<float>(varname[0]);
+    auto varVY = bpIn.InquireVariable<float>(varname[1]);
+    auto varVZ = bpIn.InquireVariable<float>(varname[2]);
+    auto varCurl = bpIn.InquireVariable<float>(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<float>(i);
+                float y = static_cast<float>(j);
+                float z = static_cast<float>(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;