Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add the CURL function to derived variables #4114

Merged
merged 31 commits into from
Apr 3, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
31 commits
Select commit Hold shift + click to select a range
759d003
Derived Expressions Curl function with correctness test
lizdulac Jan 24, 2024
20bd7c8
Fixing index bug for internal boundaries
lizdulac Jan 26, 2024
a73f535
Curl Correctness Test Cleanup
lizdulac Jan 26, 2024
dd153a4
Remove debug statements in Test
lizdulac Feb 6, 2024
1a6c011
Merging Bison3.8 parser into ADIOS (ERRORS - cannot build)
lizdulac Feb 20, 2024
98997ca
CMake Tweaks
eisenhauer Feb 20, 2024
69f54f1
whoops, missed the cmake file.
eisenhauer Feb 21, 2024
3c25ccb
Updated parser merged with adios, updated error message for missing o…
lizdulac Feb 22, 2024
f0edbf3
Merge new parser with curl implementation
lizdulac Feb 22, 2024
9cec433
Clean up
lizdulac Feb 26, 2024
9e30fb5
Removed 'var' keyword from grammar and correctness test, and minor ch…
lizdulac Feb 28, 2024
28ba751
merge with master
lizdulac Mar 27, 2024
53f783c
Set Derived_Variable flag to ON (for CI testing)
lizdulac Mar 27, 2024
984d344
Remove unnecessary files
lizdulac Mar 28, 2024
f7d3828
Cleanup curl code
lizdulac Mar 28, 2024
3360e38
Merge branch 'master' into curl
lizdulac Mar 28, 2024
73f2dde
Formatting
lizdulac Mar 28, 2024
aac2ec2
Formatting
lizdulac Mar 28, 2024
7061011
Restore examples CMakeLists.txt file
lizdulac Mar 28, 2024
d6c4a57
Restore cmake compile flags
lizdulac Mar 28, 2024
b69f7fb
Unused variable curlV in testing
lizdulac Mar 28, 2024
e8b4b8a
Implicit type conversion warnings
lizdulac Apr 1, 2024
bd60f50
Formatting
lizdulac Apr 1, 2024
488bfbd
Use float specific math, remove infinity check
lizdulac Apr 1, 2024
5162c52
Merge branch 'master' into curl
lizdulac Apr 1, 2024
0156489
Convert size_t iterator to float math
lizdulac Apr 1, 2024
bb5bd58
Default to DerivedVariables OFF
lizdulac Apr 1, 2024
b7a44c9
Default to DerivedVariables ON
lizdulac Apr 1, 2024
f136d21
Merge branch 'master' into curl
lizdulac Apr 2, 2024
7d820fb
Change the way to check for dimension equality and function rename
anagainaru Apr 3, 2024
e9699a3
format
anagainaru Apr 3, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just checking whether we are turning ON derived variables for 2.10.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We are not, I am reviewing this today and we turn it off before we merge

adios_option(PIP "Enable support for pip packaging" OFF)

option(ADIOS2_Blosc2_PREFER_SHARED "Prefer shared Blosc2 libraries" ON)
Expand Down
5 changes: 4 additions & 1 deletion source/adios2/toolkit/derived/ExprHelper.h
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,8 @@ enum ExpressionOperator
OP_ADD,
OP_SQRT,
OP_POW,
OP_MAGN
OP_MAGN,
OP_CURL
};

struct OperatorProperty
Expand All @@ -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 = {
Expand All @@ -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; }
Expand Down
108 changes: 106 additions & 2 deletions source/adios2/toolkit/derived/Function.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -54,15 +54,119 @@ Dims SameDimsFunc(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());
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<std::invalid_argument>("Derived", "Function", "SameDimFunc",
helper::Throw<std::invalid_argument>("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<Dims> 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<std::invalid_argument>("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 <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)>);

Expand Down
6 changes: 6 additions & 0 deletions source/adios2/toolkit/derived/Function.h
Original file line number Diff line number Diff line change
Expand Up @@ -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 CurlDimsFunc(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, CurlDimsFunc}},
{adios2::detail::ExpressionOperator::OP_MAGN, {MagnitudeFunc, SameDimsFunc}}};

template <class T>
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -1903,4 +1903,3 @@ adios2::detail::ASTDriver::parse (const std::string input)
parse.set_debug_level (trace_parsing);
parse ();
}

152 changes: 151 additions & 1 deletion testing/adios2/derived/TestBPDerivedCorrectness.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,6 @@

#include <cmath>
#include <iostream>
#include <limits>
#include <numeric>
#include <random>
#include <stdexcept>
Expand Down Expand Up @@ -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;
Expand Down
Loading