Skip to content

Commit

Permalink
Extend Expressions: 3D Cross product (ornladios#4319)
Browse files Browse the repository at this point in the history 
* Create cross product

* switch cross product dims function to new format

* formatting

* formatting

* remove unused variable
  • Loading branch information
@lizdulac
lizdulac authored Aug 28, 2024
1 parent b212c62 commit 5020069
Show file tree
Hide file tree
Showing 5 changed files with 156 additions and 34 deletions.
38 changes: 20 additions & 18 deletions source/adios2/toolkit/derived/Expression.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,27 +36,28 @@ const std::map<ExpressionOperator, OperatorProperty> op_property = {
{ExpressionOperator::OP_ASIN, {"ASIN", false}},
{ExpressionOperator::OP_ACOS, {"ACOS", false}},
{ExpressionOperator::OP_ATAN, {"ATAN", false}},
{ExpressionOperator::OP_CURL, {"CURL", false}},
{ExpressionOperator::OP_MAGN, {"MAGNITUDE", false}}};
{ExpressionOperator::OP_MAGN, {"MAGNITUDE", false}},
{ExpressionOperator::OP_CROSS, {"CROSS", false}},
{ExpressionOperator::OP_CURL, {"CURL", false}}};

const std::map<std::string, ExpressionOperator> string_to_op = {
{"ALIAS", ExpressionOperator::OP_ALIAS}, /* Parser-use only */
{"PATH", ExpressionOperator::OP_PATH}, /* Parser-use only */
{"NUM", ExpressionOperator::OP_NUM}, /* Parser-use only */
{"INDEX", ExpressionOperator::OP_INDEX}, {"+", ExpressionOperator::OP_ADD},
{"add", ExpressionOperator::OP_ADD}, {"ADD", ExpressionOperator::OP_ADD},
{"-", ExpressionOperator::OP_SUBTRACT}, {"SUBTRACT", ExpressionOperator::OP_SUBTRACT},
{"/", ExpressionOperator::OP_DIV}, {"divide", ExpressionOperator::OP_DIV},
{"DIVIDE", ExpressionOperator::OP_DIV}, {"*", ExpressionOperator::OP_MULT},
{"multiply", ExpressionOperator::OP_MULT}, {"MULTIPLY", ExpressionOperator::OP_MULT},
{"SQRT", ExpressionOperator::OP_SQRT}, {"sqrt", ExpressionOperator::OP_SQRT},
{"pow", ExpressionOperator::OP_POW}, {"POW", ExpressionOperator::OP_POW},
{"sin", ExpressionOperator::OP_SIN}, {"cos", ExpressionOperator::OP_COS},
{"tan", ExpressionOperator::OP_TAN}, {"asin", ExpressionOperator::OP_ASIN},
{"acos", ExpressionOperator::OP_ACOS}, {"atan", ExpressionOperator::OP_ATAN},
{"^", ExpressionOperator::OP_POW}, {"CURL", ExpressionOperator::OP_CURL},
{"curl", ExpressionOperator::OP_CURL}, {"MAGNITUDE", ExpressionOperator::OP_MAGN},
{"magnitude", ExpressionOperator::OP_MAGN}};
{"INDEX", ExpressionOperator::OP_INDEX}, {"+", ExpressionOperator::OP_ADD},
{"add", ExpressionOperator::OP_ADD}, {"ADD", ExpressionOperator::OP_ADD},
{"-", ExpressionOperator::OP_SUBTRACT}, {"SUBTRACT", ExpressionOperator::OP_SUBTRACT},
{"/", ExpressionOperator::OP_DIV}, {"divide", ExpressionOperator::OP_DIV},
{"DIVIDE", ExpressionOperator::OP_DIV}, {"*", ExpressionOperator::OP_MULT},
{"multiply", ExpressionOperator::OP_MULT}, {"MULTIPLY", ExpressionOperator::OP_MULT},
{"SQRT", ExpressionOperator::OP_SQRT}, {"sqrt", ExpressionOperator::OP_SQRT},
{"pow", ExpressionOperator::OP_POW}, {"POW", ExpressionOperator::OP_POW},
{"sin", ExpressionOperator::OP_SIN}, {"cos", ExpressionOperator::OP_COS},
{"tan", ExpressionOperator::OP_TAN}, {"asin", ExpressionOperator::OP_ASIN},
{"acos", ExpressionOperator::OP_ACOS}, {"atan", ExpressionOperator::OP_ATAN},
{"^", ExpressionOperator::OP_POW}, {"magnitude", ExpressionOperator::OP_MAGN},
{"MAGNITUDE", ExpressionOperator::OP_MAGN}, {"cross", ExpressionOperator::OP_CROSS},
{"curl", ExpressionOperator::OP_CURL}, {"CURL", ExpressionOperator::OP_CURL}};

inline std::string get_op_name(ExpressionOperator op) { return op_property.at(op).name; }

Expand Down Expand Up @@ -151,9 +152,10 @@ std::map<adios2::detail::ExpressionOperator, OperatorFunctions> OpFunctions = {
{adios2::detail::ExpressionOperator::OP_ASIN, {AsinFunc, SameDimsFunc, FloatTypeFunc}},
{adios2::detail::ExpressionOperator::OP_ACOS, {AcosFunc, SameDimsFunc, FloatTypeFunc}},
{adios2::detail::ExpressionOperator::OP_ATAN, {AtanFunc, SameDimsFunc, FloatTypeFunc}},
{adios2::detail::ExpressionOperator::OP_CURL, {Curl3DFunc, CurlDimsFunc, SameTypeFunc}},
{adios2::detail::ExpressionOperator::OP_MAGN,
{MagnitudeFunc, SameDimsWithAgrFunc, SameTypeFunc}}};
{MagnitudeFunc, SameDimsWithAgrFunc, SameTypeFunc}},
{adios2::detail::ExpressionOperator::OP_CROSS, {Cross3DFunc, Cross3DDimsFunc, SameTypeFunc}},
{adios2::detail::ExpressionOperator::OP_CURL, {Curl3DFunc, CurlDimsFunc, SameTypeFunc}}};

Expression::Expression(std::string string_exp)
: m_Shape({0}), m_Start({0}), m_Count({0}), ExprString(string_exp)
Expand Down
1 change: 1 addition & 0 deletions source/adios2/toolkit/derived/Expression.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ enum ExpressionOperator
OP_ACOS,
OP_ATAN,
OP_MAGN,
OP_CROSS,
OP_CURL
};
}
Expand Down
70 changes: 70 additions & 0 deletions source/adios2/toolkit/derived/Function.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,6 +66,20 @@ inline size_t returnIndex(size_t x, size_t y, size_t z, const size_t dims[3])
return z + y * dims[2] + x * dims[2] * dims[1];
}

template <class T>
T *ApplyCross3D(const T *Ax, const T *Ay, const T *Az, const T *Bx, const T *By, const T *Bz,
const size_t dataSize)
{
T *data = (T *)malloc(dataSize * sizeof(T) * 3);
for (size_t i = 0; i < dataSize; ++i)
{
data[3 * i] = (Ay[i] * Bz[i]) - (Az[i] * By[i]);
data[3 * i + 1] = (Ax[i] * Bz[i]) - (Az[i] * Bx[i]);
data[3 * i + 2] = (Ax[i] * By[i]) - (Ay[i] * Bx[i]);
}
return data;
}

template <class T>
T *ApplyCurl(const T *input1, const T *input2, const T *input3, const size_t dims[3])
{
Expand Down Expand Up @@ -586,6 +600,41 @@ DerivedData MagnitudeFunc(std::vector<DerivedData> inputData, DataType type)
return DerivedData();
}

DerivedData Cross3DFunc(const std::vector<DerivedData> inputData, DataType type)
{
PERFSTUBS_SCOPED_TIMER("derived::Function::Cross3DFunc");
if (inputData.size() != 6)
{
helper::Throw<std::invalid_argument>("Derived", "Function", "Cross3DFunc",
"Invalid number of arguments passed to Cross3DFunc");
}
size_t dataSize = std::accumulate(std::begin(inputData[0].Count), std::end(inputData[0].Count),
1, std::multiplies<size_t>());

DerivedData cross;
cross.Start = inputData[0].Start;
cross.Start.push_back(0);
cross.Count = inputData[0].Count;
cross.Count.push_back(3);

#define declare_type_cross(T) \
if (type == helper::GetDataType<T>()) \
{ \
T *Ax = (T *)inputData[0].Data; \
T *Ay = (T *)inputData[1].Data; \
T *Az = (T *)inputData[2].Data; \
T *Bx = (T *)inputData[3].Data; \
T *By = (T *)inputData[4].Data; \
T *Bz = (T *)inputData[5].Data; \
cross.Data = detail::ApplyCross3D(Ax, Ay, Az, Bx, By, Bz, dataSize); \
return cross; \
}
ADIOS2_FOREACH_ATTRIBUTE_PRIMITIVE_STDTYPE_1ARG(declare_type_cross)
helper::Throw<std::invalid_argument>("Derived", "Function", "Cross3DFunc",
"Invalid variable types");
return DerivedData();
}

/*
* Input: 3D vector field F(x,y,z)= {F1(x,y,z), F2(x,y,z), F3(x,y,z)}
*
Expand Down Expand Up @@ -660,6 +709,27 @@ std::tuple<Dims, Dims, Dims> SameDimsWithAgrFunc(std::vector<std::tuple<Dims, Di
return {outStart, outCount, outShape};
}

std::tuple<Dims, Dims, Dims> Cross3DDimsFunc(std::vector<std::tuple<Dims, Dims, Dims>> input)
{
// check that all dimenstions are the same
if (input.size() > 1)
{
auto first_element = input[0];
bool dim_are_equal = std::all_of(
input.begin() + 1, input.end(),
[&first_element](std::tuple<Dims, Dims, Dims> x) { return x == first_element; });
if (!dim_are_equal)
helper::Throw<std::invalid_argument>("Derived", "Function", "Cross3DDimsFunc",
"Invalid variable dimensions");
}
// return original dimensions with added dimension of number of inputs
std::tuple<Dims, Dims, Dims> output = input[0];
std::get<0>(output).push_back(0);
std::get<1>(output).push_back(3);
std::get<2>(output).push_back(3);
return output;
}

// Input Dims are the same, output is combination of all inputs
std::tuple<Dims, Dims, Dims> CurlDimsFunc(std::vector<std::tuple<Dims, Dims, Dims>> input)
{
Expand Down
2 changes: 2 additions & 0 deletions source/adios2/toolkit/derived/Function.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,10 +20,12 @@ DerivedData DivFunc(std::vector<DerivedData> input, DataType type);
DerivedData SqrtFunc(std::vector<DerivedData> input, DataType type);
DerivedData PowFunc(std::vector<DerivedData> input, DataType type);
DerivedData MagnitudeFunc(std::vector<DerivedData> input, DataType type);
DerivedData Cross3DFunc(std::vector<DerivedData> input, DataType type);
DerivedData Curl3DFunc(std::vector<DerivedData> input, DataType type);

std::tuple<Dims, Dims, Dims> SameDimsFunc(std::vector<std::tuple<Dims, Dims, Dims>> input);
std::tuple<Dims, Dims, Dims> SameDimsWithAgrFunc(std::vector<std::tuple<Dims, Dims, Dims>> input);
std::tuple<Dims, Dims, Dims> Cross3DDimsFunc(std::vector<std::tuple<Dims, Dims, Dims>> input);
std::tuple<Dims, Dims, Dims> CurlDimsFunc(std::vector<std::tuple<Dims, Dims, Dims>> input);

DataType SameTypeFunc(DataType input);
Expand Down
79 changes: 63 additions & 16 deletions testing/adios2/derived/TestBPDerivedCorrectness.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -320,7 +320,7 @@ TEST_P(DerivedCorrectnessP, TrigCorrectnessTest)
bpFileReader.Close();
}

TEST_P(DerivedCorrectnessP, MagCorrectnessTest)
TEST_P(DerivedCorrectnessP, VectorCorrectnessTest)
{
const size_t Nx = 2, Ny = 3, Nz = 10;
const size_t steps = 2;
Expand All @@ -334,30 +334,50 @@ TEST_P(DerivedCorrectnessP, MagCorrectnessTest)
std::vector<float> simArray1(Nx * Ny * Nz);
std::vector<float> simArray2(Nx * Ny * Nz);
std::vector<float> simArray3(Nx * Ny * Nz);
std::vector<float> simArray4(Nx * Ny * Nz);
std::vector<float> simArray5(Nx * Ny * Nz);
std::vector<float> simArray6(Nx * Ny * Nz);
for (size_t i = 0; i < Nx * Ny * Nz; ++i)
{
simArray1[i] = distribution(generator);
simArray2[i] = distribution(generator);
simArray3[i] = distribution(generator);
simArray4[i] = distribution(generator);
simArray5[i] = distribution(generator);
simArray6[i] = distribution(generator);
}

adios2::ADIOS adios;
adios2::IO bpOut = adios.DeclareIO("BPWriteExpression");
std::vector<std::string> varname = {"sim2/Ux", "sim2/Uy", "sim2/Uz"};
std::string derivedname = "derived/magU";
std::vector<std::string> varname = {"sim2/Ux", "sim2/Uy", "sim2/Uz",
"sim2/Vx", "sim2/Vy", "sim2/Vz"};
const std::string derMagName = "derived/magU";
const std::string derCrossName = "derived/crossU";

auto Ux = bpOut.DefineVariable<float>(varname[0], {Nx, Ny, Nz}, {0, 0, 0}, {Nx, Ny, Nz});
auto Uy = bpOut.DefineVariable<float>(varname[1], {Nx, Ny, Nz}, {0, 0, 0}, {Nx, Ny, Nz});
auto Uz = bpOut.DefineVariable<float>(varname[2], {Nx, Ny, Nz}, {0, 0, 0}, {Nx, Ny, Nz});
auto Vx = bpOut.DefineVariable<float>(varname[3], {Nx, Ny, Nz}, {0, 0, 0}, {Nx, Ny, Nz});
auto Vy = bpOut.DefineVariable<float>(varname[4], {Nx, Ny, Nz}, {0, 0, 0}, {Nx, Ny, Nz});
auto Vz = bpOut.DefineVariable<float>(varname[5], {Nx, Ny, Nz}, {0, 0, 0}, {Nx, Ny, Nz});
// clang-format off
bpOut.DefineDerivedVariable(derivedname,
bpOut.DefineDerivedVariable(derMagName,
"x =" + varname[0] + " \n"
"y =" + varname[1] + " \n"
"z =" + varname[2] + " \n"
"magnitude(x,y,z)",
mode);
bpOut.DefineDerivedVariable(derCrossName,
"Ux =" + varname[0] + " \n"
"Uy =" + varname[1] + " \n"
"Uz =" + varname[2] + " \n"
"Vx =" + varname[3] + " \n"
"Vy =" + varname[4] + " \n"
"Vz =" + varname[5] + " \n"
"cross(Ux, Uy, Uz, Vx, Vy, Vz)",
mode);
// clang-format on
std::string filename = "expMagnitude.bp";
std::string filename = "expVector.bp";
adios2::Engine bpFileWriter = bpOut.Open(filename, adios2::Mode::Write);

for (size_t i = 0; i < steps; i++)
Expand All @@ -366,6 +386,9 @@ TEST_P(DerivedCorrectnessP, MagCorrectnessTest)
bpFileWriter.Put(Ux, simArray1.data());
bpFileWriter.Put(Uy, simArray2.data());
bpFileWriter.Put(Uz, simArray3.data());
bpFileWriter.Put(Vx, simArray4.data());
bpFileWriter.Put(Vy, simArray5.data());
bpFileWriter.Put(Vz, simArray6.data());
bpFileWriter.EndStep();
}
bpFileWriter.Close();
Expand All @@ -376,28 +399,52 @@ TEST_P(DerivedCorrectnessP, MagCorrectnessTest)
std::vector<float> readUx;
std::vector<float> readUy;
std::vector<float> readUz;
std::vector<float> readVx;
std::vector<float> readVy;
std::vector<float> readVz;
std::vector<float> readMag;
std::vector<float> readCross;

float calcM;
float calcDerived;
float epsilon = (float)0.01;
for (size_t i = 0; i < steps; i++)
{
bpFileReader.BeginStep();
auto varx = bpIn.InquireVariable<float>(varname[0]);
auto vary = bpIn.InquireVariable<float>(varname[1]);
auto varz = bpIn.InquireVariable<float>(varname[2]);
auto varmag = bpIn.InquireVariable<float>(derivedname);

bpFileReader.Get(varx, readUx);
bpFileReader.Get(vary, readUy);
bpFileReader.Get(varz, readUz);
auto varUx = bpIn.InquireVariable<float>(varname[0]);
auto varUy = bpIn.InquireVariable<float>(varname[1]);
auto varUz = bpIn.InquireVariable<float>(varname[2]);
auto varVx = bpIn.InquireVariable<float>(varname[3]);
auto varVy = bpIn.InquireVariable<float>(varname[4]);
auto varVz = bpIn.InquireVariable<float>(varname[5]);
auto varmag = bpIn.InquireVariable<float>(derMagName);
auto varcross = bpIn.InquireVariable<float>(derCrossName);

bpFileReader.Get(varUx, readUx);
bpFileReader.Get(varUy, readUy);
bpFileReader.Get(varUz, readUz);
bpFileReader.Get(varVx, readVx);
bpFileReader.Get(varVy, readVy);
bpFileReader.Get(varVz, readVz);
bpFileReader.Get(varmag, readMag);
bpFileReader.Get(varcross, readCross);
bpFileReader.EndStep();

// Magnitude
for (size_t ind = 0; ind < Nx * Ny * Nz; ++ind)
{
calcDerived =
(float)sqrt(pow(readUx[ind], 2) + pow(readUy[ind], 2) + pow(readUz[ind], 2));
EXPECT_TRUE(fabs(calcDerived - readMag[ind]) < epsilon);
}
// Cross Product
for (size_t ind = 0; ind < Nx * Ny * Nz; ++ind)
{
calcM = (float)sqrt(pow(readUx[ind], 2) + pow(readUy[ind], 2) + pow(readUz[ind], 2));
EXPECT_TRUE(fabs(calcM - readMag[ind]) < epsilon);
calcDerived = (readUy[ind] * readVz[ind]) - (readUz[ind] * readVy[ind]);
EXPECT_TRUE(fabs(calcDerived - readCross[3 * ind]) < epsilon);
calcDerived = (readUx[ind] * readVz[ind]) - (readUz[ind] * readVx[ind]);
EXPECT_TRUE(fabs(calcDerived - readCross[3 * ind + 1]) < epsilon);
calcDerived = (readUx[ind] * readVy[ind]) - (readUy[ind] * readVx[ind]);
EXPECT_TRUE(fabs(calcDerived - readCross[3 * ind + 2]) < epsilon);
}
}
bpFileReader.Close();
Expand Down

0 comments on commit 5020069

Please sign in to comment.