OilfieldToolsNet/IFEM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

added: (Tensor|Vec)Func::timeDerivative

Browse Source 
returns the time derivative of the function
This commit is contained in:
Arne Morten Kvarving 2023-09-21 13:36:44 +02:00
parent dbcbede16c
commit c22de7be5a
6 changed files with 190 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
src/Utility/ExprFunctions.C
View File

@ -499,6 +499,19 @@ EvalMultiFunction<ParentFunc, Ret>::evalGradient (const Vec3& X) const
}
template <class ParentFunc, class Ret>
std::vector<Real>
EvalMultiFunction<ParentFunc, Ret>::evalTimeDerivative (const Vec3& X) const
{
std::vector<Real> result;
result.reserve(this->ncmp);
for (const std::unique_ptr<EvalFunction>& f : this->p)
result.push_back(f->deriv(X,4));
return result;
}
template class EvalMultiFunction<VecFunc,Vec3>;
template class EvalMultiFunction<TensorFunc,Tensor>;
template class EvalMultiFunction<STensorFunc,SymmTensor>;

3
src/Utility/ExprFunctions.h
View File

@ -226,6 +226,9 @@ protected:
//! \brief Returns the gradient of the function as a 1D array.
std::vector<Real> evalGradient(const Vec3& X) const override;
//! \brief Returns the time derivatives of the function as a 1D array.
std::vector<Real> evalTimeDerivative(const Vec3& X) const override;
};
//! Vector-valued function expression

10
src/Utility/Function.h
View File

@ -112,6 +112,8 @@ namespace utl
protected:
//! \brief Returns the gradient of the function as a 1D array.
virtual std::vector<Real> evalGradient(const Vec3&) const { return {}; }
//! \brief Returns the time derivatives of the function as a 1D array.
virtual std::vector<Real> evalTimeDerivative(const Vec3&) const { return {}; }
Result zero; //!< Return value for default implementations of derivatives
};
@ -220,6 +222,9 @@ public:
//! \brief Returns a representative scalar equivalent of the function value.
virtual Real getScalarValue(const Vec3& X) const { return this->evaluate(X); }
//! \brief Returns the time derivative of the function.
Real timeDerivative(const Vec3& X) const { return this->deriv(X,4); }
};
@ -262,6 +267,9 @@ public:
result = this->evalGradient(X);
return result;
}
//! \brief Evaluates time derivatives of the function.
Vec3 timeDerivative(const Vec3& X) const { return this->evalTimeDerivative(X); }
};
@ -276,7 +284,7 @@ public:
virtual bool isNormalPressure() const { return false; }
//! \brief Returns the time-derivative of the function.
virtual Vec3 deriv(const Vec3&, const Vec3&) const { return Vec3(); }
virtual Vec3 timeDerivative(const Vec3&, const Vec3&) const { return Vec3(); }
};
#endif

18
src/Utility/TensorFunction.C
View File

@ -23,6 +23,15 @@ utl::matrix3d<Real> TensorFunc::gradient(const Vec3& X) const
}
Tensor TensorFunc::timeDerivative(const Vec3& X) const
{
const size_t nsd = sqrt(ncmp);
Tensor result(nsd);
result = this->evalTimeDerivative(X);
return result;
}
size_t STensorFunc::index(size_t nsd, size_t i, size_t j) const
{
if (i == j)
@ -50,3 +59,12 @@ utl::matrix3d<Real> STensorFunc::gradient(const Vec3& X) const
return result;
}
SymmTensor STensorFunc::timeDerivative(const Vec3& X) const
{
const size_t nsd = ncmp > 5 ? 3 : (ncmp > 2 ? 2 : 1);
SymmTensor result(nsd, ncmp == 4);
result = this->evalTimeDerivative(X);
return result;
}

6
src/Utility/TensorFunction.h
View File

@ -53,6 +53,9 @@ public:
//! \brief Evaluates first derivatives of the function.
utl::matrix3d<Real> gradient(const Vec3& X) const;
//! \brief Evaluates time derivatives of the function.
Tensor timeDerivative(const Vec3& X) const;
};
@ -95,6 +98,9 @@ public:
//! \brief Evaluates first derivatives of the function.
utl::matrix3d<Real> gradient(const Vec3& X) const;
//! \brief Evaluates time derivatives of the function.
SymmTensor timeDerivative(const Vec3& X) const;
};
#endif

142
src/Utility/Test/TestFunctions.C
View File

@ -320,6 +320,61 @@ TEST(TestVecFuncExpr, NumDimensions)
}
TEST(TestVecFuncExpr, TimeDerivative)
{
const char* g = "sin(x)*sin(y)*sin(z)*sin(t) | x*x*y*y*z*t*t | exp(-2*t)";
const char* g_t = "sin(x)*sin(y)*sin(z)*cos(t) | x*x*y*y*z*2*t | -2*exp(-2*t)";
VecFuncExpr f(g);
f.addDerivative(g_t,"",4);
EXPECT_FALSE(f.isConstant());
for (double t : {0.1, 0.2, 0.3})
for (double x : {0.1, 0.2, 0.3})
for (double y : {0.5, 0.6, 0.7})
for (double z : {0.8, 0.9, 1.0}) {
const Vec4 X(x,y,z,t);
const Vec3 r(sin(x)*sin(y)*sin(z)*cos(t),
x*x*y*y*z*2*t,
-2*exp(-2*t));
const Vec3 dt = f.timeDerivative(X);
EXPECT_DOUBLE_EQ(dt[0], r[0]);
EXPECT_DOUBLE_EQ(dt[1], r[1]);
EXPECT_DOUBLE_EQ(dt[2], r[2]);
}
}
TEST(TestVecFuncExpr, TimeDerivativeFD)
{
const char* g = "sin(x)*sin(y)*sin(z)*sin(t) | x*x*y*y*z*t*t | exp(-2*t)";
const double eps = 1e-6;
VecFuncExpr f(g,"",1e-8,eps);
EXPECT_FALSE(f.isConstant());
for (double t : {0.1, 0.2, 0.3})
for (double x : {0.1, 0.2, 0.3})
for (double y : {0.5, 0.6, 0.7})
for (double z : {0.8, 0.9, 1.0}) {
const Vec4 X(x,y,z,t);
const Vec4 Xp(x,y,z,t + 0.5*eps);
const Vec4 Xm(x,y,z,t - 0.5*eps);
Vec3 r(sin(x)*sin(y)*sin(z)*(sin(Xp.t) - sin(Xm.t)),
x*x*y*y*z*(Xp.t*Xp.t - Xm.t*Xm.t),
exp(-2*Xp.t) - exp(-2*Xm.t));
r *= 1.0 / eps;
const Vec3 dt = f.timeDerivative(X);
EXPECT_NEAR(dt[0], r[0], 1e-8);
EXPECT_NEAR(dt[1], r[1], 1e-8);
EXPECT_NEAR(dt[2], r[2], 1e-8);
}
}
TEST(TestTensorFunc, Evaluate)
{
const char* func = "sin(x) | cos (y) | exp(z) | sin(x)*cos(y)";
@ -515,7 +570,7 @@ TEST(TestTensorFunction, Gradient3DFD)
EXPECT_NEAR(grad(i,j,d), r(i,j,d), 1e-8);
}
}
}
}
}
@ -647,6 +702,39 @@ TEST(TestTensorFunction, Hessian3D)
}
TEST(TestTensorFunction, TimeDerivative)
{
const char* g = "sin(x)*sin(y)*sin(z)*sin(t) | x*x*y*y*z*t*t | exp(x)*exp(2*y)*z*z*t |"
"exp(-2*x)*exp(y)*z*sin(t) | x*y*z*t*t | x*y*z*z*t |"
"x*sin(t) | y*t*t | z*t";
const char* g_t = "sin(x)*sin(y)*sin(z)*cos(t) | x*x*y*y*z*2*t | exp(x)*exp(2*y)*z*z |"
"exp(-2*x)*exp(y)*z*cos(t) | x*y*z*2*t | x*y*z*z |"
"x*cos(t) | y*2*t | z";
TensorFuncExpr f(g);
f.addDerivative(g_t,"",4);
EXPECT_FALSE(f.isConstant());
for (double t : {0.1, 0.2, 0.3})
for (double x : {0.1, 0.2, 0.3})
for (double y : {0.5, 0.6, 0.7})
for (double z : {0.8, 0.9, 1.0}) {
const Vec4 X(x,y,z,t);
const Tensor r({sin(x)*sin(y)*sin(z)*cos(t), x*x*y*y*z*2*t, exp(x)*exp(2*y)*z*z,
exp(-2*x)*exp(y)*z*cos(t), x*y*z*2*t, x*y*z*z,
x*cos(t), y*2*t, z});
const Tensor dt = f.timeDerivative(X);
for (size_t i = 1; i <= 3; ++i)
for (size_t j = 1; j <= 3; ++j)
EXPECT_DOUBLE_EQ(dt(i,j), r(i,j));
}
}
TEST(TestTensorFuncExpr, NumDimensions)
{
const char* func1 = "x";
@ -953,6 +1041,58 @@ TEST(TestSTensorFuncExpr, NumDimensions)
}
TEST(TestSTensorFunction, TimeDerivative)
{
const char* g = "sin(x)*sin(y)*sin(t) | exp(x)*exp(2*y)*exp(-4*t) | x*x*y*y*t*t";
const char* g_t = "sin(x)*sin(y)*cos(t) | exp(x)*exp(2*y)*-4*exp(-4*t) | x*x*y*y*2*t";
STensorFuncExpr f(g);
f.addDerivative(g_t,"",4);
EXPECT_FALSE(f.isConstant());
for (double t : {0.1, 0.2, 0.3})
for (double x : {0.1, 0.2, 0.3})
for (double y : {0.5, 0.6, 0.7}) {
const Vec4 X(x,y,0,t);
const SymmTensor r({sin(x)*sin(y)*cos(t), exp(x)*exp(2*y)*-4*exp(-4*t), x*x*y*y*2*t});
const SymmTensor dt = f.timeDerivative(X);
for (size_t i = 1; i <= 2; ++i)
for (size_t j = 1; j <= 2; ++j)
EXPECT_DOUBLE_EQ(dt(i,j), r(i,j));
}
}
TEST(TestSTensorFunction, TimeDerivativeFD)
{
const char* g = "sin(x)*sin(y)*sin(t) | exp(x)*exp(2*y)*exp(-4*t) | x*x*y*y*t*t";
const double eps = 1e-6;
STensorFuncExpr f(g,"",1e-8,eps);
EXPECT_FALSE(f.isConstant());
for (double t : {0.1, 0.2, 0.3})
for (double x : {0.1, 0.2, 0.3})
for (double y : {0.5, 0.6, 0.7}) {
const Vec4 X(x,y,0,t);
const Vec4 Xp(x,y,0,t + 0.5*eps);
const Vec4 Xm(x,y,0,t - 0.5*eps);
SymmTensor r({sin(x)*sin(y)*(sin(Xp.t) - sin(Xm.t)),
exp(x)*exp(2*y)*(exp(-4*Xp.t) - exp(-4*Xm.t)),
x*x*y*y*(Xp.t*Xp.t - Xm.t*Xm.t)});
r *= 1.0 / eps;
const SymmTensor dt = f.timeDerivative(X);
for (size_t i = 1; i <= 2; ++i)
for (size_t j = 1; j <= 2; ++j)
EXPECT_NEAR(dt(i,j), r(i,j), 1e-8);
}
}
TEST(TestEvalFunction, ExtraParam)
{
EvalFunction f("x*foo");