Magnitude of Curl

CMakeLists.txt

@@ -186,7 +186,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)
 adios_option(XRootD     "Enable support for XRootD" AUTO)
 adios_option(KVCACHE    "Enable support for KVCache" AUTO)

examples/CMakeLists.txt

@@ -11,3 +11,7 @@ add_subdirectory(useCases)
 if(ADIOS2_USE_Campaign)
     add_subdirectory(campaign)
 endif()
+if(ADIOS2_USE_Derived_Variable)
+    add_subdirectory(derived)
+endif()
+

examples/derived/CMakeLists.txt

@@ -0,0 +1,9 @@
+#------------------------------------------------------------------------------#
+# Distributed under the OSI-approved Apache License, Version 2.0.  See
+# accompanying file Copyright.txt for details.
+#------------------------------------------------------------------------------#
+
+add_executable(adios2_read_write_derived_serial read_write_derived.cpp)
+target_link_libraries(adios2_read_write_derived_serial adios2::cxx11 adios2_core)
+add_executable(adios2_write_derived_serial write_derived.cpp)
+target_link_libraries(adios2_write_derived_serial adios2::cxx11 adios2_core)

examples/derived/read_write_derived.cpp

@@ -0,0 +1,74 @@
+#include <cstdint>
+#include <cstring>
+
+#include <cmath>
+#include <iostream>
+#include <numeric>
+#include <random>
+#include <stdexcept>
+#include <vector>
+
+#include <adios2.h>
+
+int main(int argc, char **argv)
+{
+    const size_t Nx = 25, Ny = 70, Nz = 13;
+
+    // Application variables
+    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);
+            }
+        }
+    }
+
+    adios2::ADIOS adios;
+    adios2::IO bpOut = adios.DeclareIO("BPWriteExpression");
+
+    adios2::IO bpIn = adios.DeclareIO("BPReadCurlExpression");
+    std::string filename = "expCurl.bp";
+    std::string derivedname = "derived/curlV";
+    adios2::Engine bpFileReader = bpIn.Open(filename, adios2::Mode::Read);
+
+    std::vector<float> readCurl;
+
+    bpFileReader.BeginStep();
+    auto varCurl = bpIn.InquireVariable<float>(derivedname);
+    bpFileReader.Get(varCurl, readCurl);
+    bpFileReader.EndStep();
+
+    auto curlV =
+        bpOut.DefineVariable<float>("copied/curlV", {Nx, Ny, Nz, 3}, {0, 0, 0, 0}, {Nx, Ny, Nz, 3});
+    // clang-format off
+    bpOut.DefineDerivedVariable("derived/magofcurl",
+                               "curlV =copied/curlV \n"
+                               "magnitude(curlV, 3)",
+                               adios2::DerivedVarType::StoreData);
+    // clang-format on
+
+    adios2::Engine bpFileWriter = bpOut.Open("expMagCopiedCurl.bp", adios2::Mode::Write);
+
+    bpFileWriter.BeginStep();
+    bpFileWriter.Put(curlV, readCurl.data());
+    bpFileWriter.EndStep();
+    bpFileWriter.Close();
+
+    bpFileReader.Close();
+
+    return 0;
+}

examples/derived/write_derived.cpp

@@ -0,0 +1,102 @@
+#include <cstdint>
+#include <cstring>
+
+#include <cmath>
+#include <iostream>
+#include <numeric>
+#include <random>
+#include <stdexcept>
+#include <vector>
+
+#include <adios2.h>
+
+int main(int argc, char **argv)
+{
+    const size_t Nn = 320;
+    const size_t Nx = Nn, Ny = Nn, Nz = Nn;
+
+    // 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");
+
+    auto VX = bpOut.DefineVariable<float>("sim/VX", {Nx, Ny, Nz}, {0, 0, 0}, {Nx, Ny, Nz});
+    auto VY = bpOut.DefineVariable<float>("sim/VY", {Nx, Ny, Nz}, {0, 0, 0}, {Nx, Ny, Nz});
+    auto VZ = bpOut.DefineVariable<float>("sim/VZ", {Nx, Ny, Nz}, {0, 0, 0}, {Nx, Ny, Nz});
+    // clang-format off
+    //*
+    bpOut.DefineDerivedVariable("derived/curlV",
+                                "Vx =sim/VX \n"
+                                "Vy =sim/VY \n"
+                                "Vz =sim/VZ \n"
+                                "curl(Vx,Vy,Vz)",
+                                adios2::DerivedVarType::StoreData);
+    //*/
+    /*
+    bpOut.DefineDerivedVariable("derived/magV",
+                                "Vx =sim/VX \n"
+                                "Vy =sim/VY \n"
+                                "Vz =sim/VZ \n"
+                                "magnitude(Vx,Vy,Vz)",
+                                adios2::DerivedVarType::StoreData);
+    */
+    /*
+    bpOut.DefineDerivedVariable("derived/addV",
+                                "Vx =sim/VX \n"
+                                "Vy =sim/VY \n"
+                                "Vx + Vy",
+                                adios2::DerivedVarType::StoreData);
+    */
+    /*
+    bpOut.DefineDerivedVariable("derived/magofcurl",
+                                "curlx =derived/curlV[::3] \n"
+                                "curly =derived/curlV[1::3] \n"
+                                "curlz =derived/curlV[2::3] \n"
+                                "magnitude(curlx,curly,curlz)",
+                                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();
+
+    std::cout << "Example complete, check " << filename << " for data" << std::endl;
+
+    return 0;
+}

source/adios2/toolkit/derived/Expression.cpp

@@ -235,7 +235,7 @@ std::string ExpressionTree::toStringExpr()
             if (!detail.indices.empty())
             {
                 result += "[ ";
-                for (std::tuple<int, int, int> idx : detail.indices)
+                for (std::tuple<size_t, size_t, size_t> idx : detail.indices)
                 {
                     result += (std::get<0>(idx) < 0 ? "" : std::to_string(std::get<0>(idx))) + ":";
                     result += (std::get<1>(idx) < 0 ? "" : std::to_string(std::get<1>(idx))) + ":";
@@ -308,9 +308,23 @@ ExpressionTree::ApplyExpression(DataType type, size_t numBlocks,
     // apply the computation operator on all blocks
     std::vector<DerivedData> outputData(numBlocks);
     auto op_fct = OpFunctions.at(detail.operation);
-    for (size_t blk = 0; blk < numBlocks; blk++)
+    // If function called over single expression with a constant,
+    // (ex: magnitude(curl(x,y,z), 3))
+    // assume user wants to extract dimension
+    if (detail.constant > 0 && sub_exprs.size() == 1)
     {
-        outputData[blk] = op_fct.ComputeFct(exprData[blk], type);
+        for (size_t blk = 0; blk < numBlocks; blk++)
+        {
+            outputData[blk] = op_fct.ComputeFct(
+                ExtractDimensionN(exprData[blk][0], type, (size_t)detail.constant), type);
+        }
+    }
+    else
+    {
+        for (size_t blk = 0; blk < numBlocks; blk++)
+        {
+            outputData[blk] = op_fct.ComputeFct(exprData[blk], type);
+        }
     }
     // deallocate intermediate data after computing the operation
     for (size_t blk = 0; blk < numBlocks; blk++)

source/adios2/toolkit/derived/Function.cpp

@@ -44,13 +44,13 @@ T *ApplyCurl(const T *input1, const T *input2, const T *input3, const size_t dim
     size_t dataSize = dims[0] * dims[1] * dims[2];
     T *data = (T *)malloc(dataSize * sizeof(float) * 3);
     size_t index = 0;
-    for (int i = 0; i < dims[0]; ++i)
+    for (int i = 0; i < (int)dims[0]; ++i)
     {
         size_t prev_i = std::max(0, i - 1), next_i = std::min((int)dims[0] - 1, i + 1);
-        for (int j = 0; j < dims[1]; ++j)
+        for (int j = 0; j < (int)dims[1]; ++j)
         {
             size_t prev_j = std::max(0, j - 1), next_j = std::min((int)dims[1] - 1, j + 1);
-            for (int k = 0; k < dims[2]; ++k)
+            for (int k = 0; k < (int)dims[2]; ++k)
             {
                 size_t prev_k = std::max(0, k - 1), next_k = std::min((int)dims[2] - 1, k + 1);
                 // curl[0] = dv3 / dy - dv2 / dz
@@ -80,21 +80,6 @@ T *ApplyCurl(const T *input1, const T *input2, const T *input3, const size_t dim
     }
     return data;
 }
-
-// types not supported for curl
-std::complex<float> *ApplyCurl(const std::complex<float> * /*input 1*/,
-                               const std::complex<float> * /*input 2*/,
-                               const std::complex<float> * /*input 3*/, const size_t[3] /*dims*/)
-{
-    return NULL;
-}
-
-std::complex<double> *ApplyCurl(const std::complex<double> * /*input 1*/,
-                                const std::complex<double> * /*input 2*/,
-                                const std::complex<double> * /*input 3*/, const size_t[3] /*dims*/)
-{
-    return NULL;
-}
 }
 
 namespace derived
@@ -112,7 +97,7 @@ DerivedData AddFunc(std::vector<DerivedData> inputData, DataType type)
                                                 [](T a, T b) { return a + b; });                   \
         return DerivedData({(void *)addValues, inputData[0].Start, inputData[0].Count});           \
     }
-    ADIOS2_FOREACH_PRIMITIVE_STDTYPE_1ARG(declare_type_add)
+    ADIOS2_FOREACH_ATTRIBUTE_PRIMITIVE_STDTYPE_1ARG(declare_type_add)
     helper::Throw<std::invalid_argument>("Derived", "Function", "AddFunc",
                                          "Invalid variable types");
     return DerivedData();
@@ -134,7 +119,7 @@ DerivedData SubtractFunc(std::vector<DerivedData> inputData, DataType type)
                 *(reinterpret_cast<T *>(inputData[0].Data) + i) - subtractValues[i];               \
         return DerivedData({(void *)subtractValues, inputData[0].Start, inputData[0].Count});      \
     }
-    ADIOS2_FOREACH_PRIMITIVE_STDTYPE_1ARG(declare_type_subtract)
+    ADIOS2_FOREACH_ATTRIBUTE_PRIMITIVE_STDTYPE_1ARG(declare_type_subtract)
     helper::Throw<std::invalid_argument>("Derived", "Function", "SubtractFunc",
                                          "Invalid variable types");
     return DerivedData();
@@ -156,7 +141,7 @@ DerivedData MagnitudeFunc(std::vector<DerivedData> inputData, DataType type)
         }                                                                                          \
         return DerivedData({(void *)magValues, inputData[0].Start, inputData[0].Count});           \
     }
-    ADIOS2_FOREACH_PRIMITIVE_STDTYPE_1ARG(declare_type_mag)
+    ADIOS2_FOREACH_ATTRIBUTE_PRIMITIVE_STDTYPE_1ARG(declare_type_mag)
     helper::Throw<std::invalid_argument>("Derived", "Function", "MagnitudeFunc",
                                          "Invalid variable types");
     return DerivedData();
@@ -196,12 +181,67 @@ DerivedData Curl3DFunc(const std::vector<DerivedData> inputData, DataType type)
         curl.Data = detail::ApplyCurl(input1, input2, input3, dims);                               \
         return curl;                                                                               \
     }
-    ADIOS2_FOREACH_PRIMITIVE_STDTYPE_1ARG(declare_type_curl)
+    ADIOS2_FOREACH_ATTRIBUTE_PRIMITIVE_STDTYPE_1ARG(declare_type_curl)
     helper::Throw<std::invalid_argument>("Derived", "Function", "Curl3DFunc",
                                          "Invalid variable types");
     return DerivedData();
 }
 
+std::vector<DerivedData> ExtractDimensionN(DerivedData inputData, DataType type, size_t dim)
+{
+    size_t num_data_sets = inputData.Count[dim];
+    size_t num_chunks = 1;
+    size_t chunk_length = 1;
+    for (size_t i = 0; i < dim; ++i)
+    {
+        num_chunks *= inputData.Count[i];
+    }
+    for (size_t i = dim + 1; i < inputData.Count.size(); ++i)
+    {
+        chunk_length *= inputData.Count[i];
+    }
+
+    Dims set_Start;
+    Dims set_Count;
+    for (size_t i = 0; i < inputData.Start.size(); ++i)
+    {
+        if (i != dim)
+        {
+            set_Start.push_back(inputData.Start[i]);
+            set_Count.push_back(inputData.Count[i]);
+        }
+    }
+
+    std::vector<DerivedData> result;
+    size_t chunk_size = chunk_length * helper::GetDataTypeSize(type);
+    // TO DO - FREE
+    for (size_t i = 0; i < num_data_sets; ++i)
+        result.push_back({malloc(num_chunks * chunk_size), set_Start, set_Count});
+
+    // How does Start factor in?
+    //  size_t data_iter = 0;
+    char *input_ptr = (char *)inputData.Data;
+    for (size_t chunk = 0; chunk < num_chunks; ++chunk)
+    {
+        for (size_t data_set = 0; data_set < num_data_sets; ++data_set)
+        {
+            char *result_ptr = (char *)(result[data_set].Data);
+            memcpy(result_ptr + (chunk * chunk_size),
+                   input_ptr + (((num_data_sets * chunk) + data_set) * chunk_size), chunk_size);
+
+            // memcpy(&(result[data_set].Data[chunk * chunk_length]), &(inputData.Data[data_iter]),
+            /*
+          for (size_t chunk_iter = 0; chunk_iter < chunk_length; ++chunk_iter)
+            {
+              result[data_set].Data[(chunk * chunk_length) + chunk_iter] =
+          inputData.Data[data_iter++];
+              }*/
+        }
+    }
+
+    return result;
+}
+
 Dims SameDimsFunc(std::vector<Dims> input)
 {
     // check that all dimenstions are the same

source/adios2/toolkit/derived/Function.h

@@ -12,6 +12,8 @@ DerivedData SubtractFunc(std::vector<DerivedData> input, DataType type);
 DerivedData MagnitudeFunc(std::vector<DerivedData> input, DataType type);
 DerivedData Curl3DFunc(std::vector<DerivedData> input, DataType type);
 
+std::vector<DerivedData> ExtractDimensionN(DerivedData inputData, DataType type, size_t dim);
+
 Dims SameDimsFunc(std::vector<Dims> input);
 Dims CurlDimsFunc(std::vector<Dims> input);
 }

source/adios2/toolkit/derived/parser/ASTDriver.cpp

@@ -84,6 +84,12 @@ void ASTDriver::createNode(std::string op_name, size_t numsubexprs)
     holding.push(node);
 }
 
+void ASTDriver::createNode(double num)
+{
+    ASTNode *node = new ASTNode("NUM", num);
+    holding.push(node);
+}
+
 void ASTDriver::createNode(std::string alias)
 {
     ASTNode *node = new ASTNode("ALIAS", alias);