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

added: allow use of automatic differentation in analytic solutions

Browse Source
This commit is contained in:
Arne Morten Kvarving 2023-09-27 09:01:20 +02:00
parent cca533fa1c
commit abb21feef7
7 changed files with 246 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
src/SIM/SIMbase.C
View File

@ -208,6 +208,9 @@ bool SIMbase::preprocess (const IntVec& ignored, bool fixDup)
static int substep = 10;
this->printHeading(substep);
if (mySol)
mySol->setupSecondarySolutions();
// Perform some sub-class specific pre-preprocessing, if any
this->preprocessA();

109
src/Utility/AnaSol.C
View File

@ -17,13 +17,51 @@
#include "Utilities.h"
#include "IFEM.h"
#include "tinyxml.h"
#include <autodiff/reverse/var.hpp>
#ifdef HAS_HDF5
#include "ProcessAdm.h"
#include "HDF5Writer.h"
#include "FieldFunctions.h"
#endif
namespace {
//! \brief Function returning the gradient of a scalar or vector function.
template<class FromFunc, class ToFunc, class Ret>
class GradientFunc : public ToFunc
{
public:
//! \brief The constructor initializes the function to use.
GradientFunc(const FromFunc& f) : func(f) {}
protected:
//! \brief Evaluates the gradient in a point.
Ret evaluate(const Vec3& X) const { return func.gradient(X); }
private:
const FromFunc& func; //!< Reference to scalar/vector function
};
//! \brief Specialization for symmetric tensors.
//! \details We need to construct a symmtensor from the returned tensor.
template<>
SymmTensor GradientFunc<VecFunc,STensorFunc,SymmTensor>::evaluate(const Vec3& X) const
{
Tensor tmp = func.gradient(X);
const size_t nsd = tmp.dim();
SymmTensor ret(nsd);
for (size_t i = 1; i <= nsd; ++i)
for (size_t j = i; j <= nsd; ++j)
ret(i,j) = tmp(i,j);
return ret;
}
}
AnaSol::AnaSol (RealFunc* s1, VecFunc* s2,
VecFunc* v1, TensorFunc* v2, STensorFunc* v3)
: vecSol(v1), vecSecSol(v2), stressSol(v3)
@ -92,8 +130,15 @@ AnaSol::AnaSol (const TiXmlElement* elem, bool scalarSol)
const char* type = elem->Attribute("type");
if (type && !strcasecmp(type,"fields"))
this->parseFieldFunctions(elem,scalarSol);
else
this->parseExpressionFunctions(elem,scalarSol);
else {
bool useAD = false;
utl::getAttribute(elem, "autodiff", useAD);
utl::getAttribute(elem, "symmetric", symmetric);
if (useAD)
this->parseExpressionFunctions<autodiff::var>(elem,scalarSol);
else
this->parseExpressionFunctions<Real>(elem,scalarSol);
}
}
@ -128,8 +173,13 @@ void AnaSol::initPatch (size_t pIdx)
}
template<class Scalar>
void AnaSol::parseExpressionFunctions (const TiXmlElement* elem, bool scalarSol)
{
using EvalF = EvalFuncSpatial<Scalar>;
using VecF = EvalMultiFunction<VecFunc,Vec3,Scalar>;
using TensorF = EvalMultiFunction<TensorFunc,Tensor,Scalar>;
using STensorF = EvalMultiFunction<STensorFunc,SymmTensor,Scalar>;
std::string variables;
const TiXmlElement* var = elem->FirstChildElement("variables");
if (var && var->FirstChild())
@ -168,16 +218,16 @@ void AnaSol::parseExpressionFunctions (const TiXmlElement* elem, bool scalarSol)
if (scalarSol)
{
if (type == "expression") {
scalSol.push_back(new EvalFunction((variables+primary).c_str()));
parseDerivatives(static_cast<EvalFunction*>(scalSol.back()),prim);
scalSol.push_back(new EvalF((variables+primary).c_str()));
parseDerivatives(static_cast<EvalF*>(scalSol.back()),prim);
} else
scalSol.push_back(utl::parseRealFunc(primary,type,false));
}
else
{
if (type == "expression") {
vecSol = new VecFuncExpr(primary,variables);
parseDerivatives(static_cast<VecFuncExpr*>(vecSol),prim);
vecSol = new VecF(primary,variables);
parseDerivatives(static_cast<VecF*>(vecSol),prim);
} else
vecSol = utl::parseVecFunc(primary, type);
}
@ -190,9 +240,9 @@ void AnaSol::parseExpressionFunctions (const TiXmlElement* elem, bool scalarSol)
utl::getAttribute(prim, "type", type);
std::string prType = (type == "expression" ? "" : "("+type+") ");
std::string primary = prim->FirstChild()->Value();
IFEM::cout <<"\tScalar Primary " << prType << "=" << primary << std::endl;
IFEM::cout <<"\tScalar Primary" << prType << "=" << primary << std::endl;
if (type == "expression")
scalSol.push_back(new EvalFunction((variables+primary).c_str()));
scalSol.push_back(new EvalF((variables+primary).c_str()));
else
scalSol.push_back(utl::parseRealFunc(primary, type, false));
prim = prim->NextSiblingElement("scalarprimary");
@ -209,18 +259,18 @@ void AnaSol::parseExpressionFunctions (const TiXmlElement* elem, bool scalarSol)
if (scalarSol)
{
if (type == "expression") {
scalSecSol.push_back(new VecFuncExpr(secondary,variables));
parseDerivatives(static_cast<VecFuncExpr*>(scalSecSol.back()),sec);
scalSecSol.push_back(new VecF(secondary,variables));
parseDerivatives(static_cast<VecF*>(scalSecSol.back()),sec);
} else
scalSecSol.push_back(utl::parseVecFunc(secondary, type));
}
else
{
if (type == "expression")
vecSecSol = new TensorFuncExpr(secondary,variables);
else
if (type == "expression") {
vecSecSol = new TensorF(secondary,variables);
parseDerivatives(static_cast<TensorF*>(vecSecSol),sec);
} else
vecSecSol = utl::parseTensorFunc(secondary, type);
parseDerivatives(static_cast<TensorFuncExpr*>(vecSecSol),sec);
}
}
@ -232,7 +282,7 @@ void AnaSol::parseExpressionFunctions (const TiXmlElement* elem, bool scalarSol)
std::string secondary = sec->FirstChild()->Value();
IFEM::cout <<"\tScalar Secondary="<< secondary << std::endl;
if (type == "expression")
scalSecSol.push_back(new VecFuncExpr(secondary,variables));
scalSecSol.push_back(new VecF(secondary,variables));
else
scalSecSol.push_back(utl::parseVecFunc(secondary, type));
sec = sec->NextSiblingElement("scalarsecondary");
@ -241,10 +291,11 @@ void AnaSol::parseExpressionFunctions (const TiXmlElement* elem, bool scalarSol)
const TiXmlElement* stress = elem->FirstChildElement("stress");
if (stress && stress->FirstChild())
{
symmetric = true;
std::string sigma = stress->FirstChild()->Value();
IFEM::cout <<"\tStress="<< sigma << std::endl;
stressSol = new STensorFuncExpr(sigma,variables);
parseDerivatives(static_cast<STensorFuncExpr*>(stressSol),stress);
stressSol = new STensorF(sigma,variables);
parseDerivatives(static_cast<STensorF*>(stressSol),stress);
}
}
@ -310,3 +361,23 @@ void AnaSol::parseFieldFunctions (const TiXmlElement* elem, bool scalarSol)
<<", no fields read."<< std::endl;
#endif
}
void AnaSol::setupSecondarySolutions ()
{
// if we are given a stress sol, no scalar secondaries should be registered.
// this is tailored to assumptions in elasticity applications.
if (!stressSol && !scalSol.empty()) {
scalSecSol.resize(scalSol.size());
for (size_t i = 0; i < scalSol.size(); ++i)
if (!scalSecSol[i])
scalSecSol[i] = new GradientFunc<RealFunc,VecFunc,Vec3>(*scalSol[i]);
}
if (vecSol) {
if (symmetric && !stressSol)
stressSol = new GradientFunc<VecFunc,STensorFunc,SymmTensor>(*vecSol);
else if (!symmetric && !vecSecSol)
vecSecSol = new GradientFunc<VecFunc,TensorFunc,Tensor>(*vecSol);
}
}

8
src/Utility/AnaSol.h
View File

@ -123,14 +123,22 @@ public:
//! \brief Sets the patch to use.
void initPatch(size_t pIdx);
//! \brief Make sure we have a secondary solution.
//! \details If none is given, we use derivation (automatic or finite difference)
//! to obtain one.
virtual void setupSecondarySolutions();
private:
//! \brief Parses expression functions from XML definition.
template<class Scalar>
void parseExpressionFunctions(const TiXmlElement* elem, bool scalarSol);
//! \brief Parses field functions from XML definition.
void parseFieldFunctions(const TiXmlElement* elem, bool scalarSol);
protected:
bool symmetric = false; //!< True to use symmetric secondary solution
std::vector<RealFunc*> scalSol; //!< Primary scalar solution fields
std::vector<VecFunc*> scalSecSol; //!< Secondary scalar solution fields

4
src/Utility/ExprFunctions.C
View File

@ -22,7 +22,7 @@
#include <autodiff/reverse/var.hpp>
template<> int EvalFunc::numError = 0;
template<> int EvalFuncScalar<Real>::numError = 0; //!< Explicit instantiation
namespace {
@ -90,7 +90,7 @@ void ExprException (const ExprEval::Exception& exc, const char* task,
std::cerr <<": Unknown exception";
}
std::cerr << std::endl;
EvalFunc::numError++;
EvalFuncScalar<Real>::numError++;
}

15
src/Utility/ExprFunctions.h
View File

@ -99,14 +99,16 @@ class EvalFuncSpatial : public RealFunc
Scalar* z; //!< Z-coordinate
Scalar* t; //!< Time
//! \brief Returns a const ref to a member.
//! \param dir One-based index to member
const Scalar& get(int dir) const
{
switch (dir) {
default:
case 1: return *x;
case 2: return *y;
case 3: return *z;
case 4: return *t;
case 1: return *x;
case 2: return *y;
case 3: return *z;
case 4: return *t;
default: return *x;
}
}
};
@ -268,4 +270,7 @@ using TensorFuncExpr = EvalMultiFunction<TensorFunc,Tensor,Real>;
//! Symmetric tensor-valued function expression
using STensorFuncExpr = EvalMultiFunction<STensorFunc,SymmTensor,Real>;
//! \brief Explicit instantiation of error flag.
template<> int EvalFunc::numError;
#endif

2
src/Utility/Functions.C
View File

@ -717,7 +717,7 @@ static const ScalarFunc* parseFunction (const char* type, char* cline, Real C)
{
cline = strtok(cline,":");
IFEM::cout << cline;
EvalFunc::numError = 0;
EvalFuncScalar<Real>::numError = 0;
ScalarFunc* sf = new EvalFunc(cline,"t",C);
if (EvalFunc::numError > 0)
{

141
src/Utility/Test/TestAnaSol.C
View File

@ -12,6 +12,8 @@
#include "AnaSol.h"
#include "Function.h"
#include "Tensor.h"
#include "TensorFunction.h"
#include "Vec3Oper.h"
#include "tinyxml.h"
@ -48,15 +50,15 @@ TEST(TestAnaSol, ParseDerivatives)
ASSERT_TRUE(v != nullptr);
Vec3 X(1.0,2.0,3.0);
EXPECT_FLOAT_EQ((*v)(X),8.0);
EXPECT_FLOAT_EQ(v->deriv(X,1),3.0);
EXPECT_FLOAT_EQ(v->deriv(X,2),4.0);
EXPECT_FLOAT_EQ(v->deriv(X,3),1.0);
EXPECT_FLOAT_EQ(v->dderiv(X,1,1),6.0);
EXPECT_FLOAT_EQ(v->dderiv(X,2,2),2.0);
EXPECT_FLOAT_EQ(v->dderiv(X,3,3),0.0);
EXPECT_FLOAT_EQ(v->dderiv(X,1,2),0.0);
EXPECT_FLOAT_EQ(v->dderiv(X,3,2),0.0);
EXPECT_DOUBLE_EQ((*v)(X),8.0);
EXPECT_DOUBLE_EQ(v->deriv(X,1),3.0);
EXPECT_DOUBLE_EQ(v->deriv(X,2),4.0);
EXPECT_DOUBLE_EQ(v->deriv(X,3),1.0);
EXPECT_DOUBLE_EQ(v->dderiv(X,1,1),6.0);
EXPECT_DOUBLE_EQ(v->dderiv(X,2,2),2.0);
EXPECT_DOUBLE_EQ(v->dderiv(X,3,3),0.0);
EXPECT_DOUBLE_EQ(v->dderiv(X,1,2),0.0);
EXPECT_DOUBLE_EQ(v->dderiv(X,3,2),0.0);
VecFunc* v2 = mySol.getScalarSecSol();
ASSERT_TRUE(v2 != nullptr);
@ -65,3 +67,124 @@ TEST(TestAnaSol, ParseDerivatives)
EXPECT_TRUE(v2->deriv(X,2) == Vec3(0.0,2.0,0.0));
EXPECT_TRUE(v2->dderiv(X,1,1) == Vec3(6.0,0.0,0.0));
}
TEST(TestAnaSol, ParseAD)
{
const char* input = "<anasol type=\"expression\" autodiff=\"true\">"
" <primary>"
" x^3+y^2+z"
" </primary>"
"</anasol>";
TiXmlDocument doc;
doc.Parse(input,nullptr,TIXML_ENCODING_UTF8);
const TiXmlElement* tag = doc.RootElement();
ASSERT_TRUE(tag != nullptr);
ASSERT_EQ(strcmp(tag->Value(),"anasol"),0);
AnaSol mySol(tag,true);
RealFunc* v = mySol.getScalarSol();
ASSERT_TRUE(v != nullptr);
mySol.setupSecondarySolutions();
Vec3 X(1.0,2.0,3.0);
const Vec3 grad1 = v->gradient(X);
const SymmTensor hess = v->hessian(X);
VecFunc* v2 = mySol.getScalarSecSol();
ASSERT_TRUE(v2 != nullptr);
const Vec3 grad2 = (*v2)(X);
EXPECT_DOUBLE_EQ((*v)(X),8.0);
EXPECT_DOUBLE_EQ(v->deriv(X,1),3.0);
EXPECT_DOUBLE_EQ(grad1[0],3.0);
EXPECT_DOUBLE_EQ(grad2[0],3.0);
EXPECT_DOUBLE_EQ(v->deriv(X,2),4.0);
EXPECT_DOUBLE_EQ(grad1[1],4.0);
EXPECT_DOUBLE_EQ(grad2[1],4.0);
EXPECT_DOUBLE_EQ(v->deriv(X,3),1.0);
EXPECT_DOUBLE_EQ(grad1[2],1.0);
EXPECT_DOUBLE_EQ(grad2[2],1.0);
EXPECT_DOUBLE_EQ(v->dderiv(X,1,1),6.0);
EXPECT_DOUBLE_EQ(hess(1,1), 6.0);
EXPECT_DOUBLE_EQ(v->dderiv(X,2,2),2.0);
EXPECT_DOUBLE_EQ(hess(2,2), 2.0);
EXPECT_DOUBLE_EQ(v->dderiv(X,3,3),0.0);
EXPECT_DOUBLE_EQ(hess(3,3), 0.0);
EXPECT_DOUBLE_EQ(v->dderiv(X,1,2),0.0);
EXPECT_DOUBLE_EQ(hess(1,2), 0.0);
EXPECT_DOUBLE_EQ(v->dderiv(X,3,2),0.0);
EXPECT_DOUBLE_EQ(hess(3,2), 0.0);
}
TEST(TestAnaSol, ParseFD)
{
const char* input = "<anasol type=\"expression\">"
" <primary>"
" x^3+y^2+z"
" </primary>"
"</anasol>";
TiXmlDocument doc;
doc.Parse(input,nullptr,TIXML_ENCODING_UTF8);
const TiXmlElement* tag = doc.RootElement();
ASSERT_TRUE(tag != nullptr);
ASSERT_EQ(strcmp(tag->Value(),"anasol"),0);
AnaSol mySol(tag,true);
RealFunc* v = mySol.getScalarSol();
ASSERT_TRUE(v != nullptr);
mySol.setupSecondarySolutions();
Vec3 X(1.0,2.0,3.0);
const Vec3 grad1 = v->gradient(X);
VecFunc* v2 = mySol.getScalarSecSol();
ASSERT_TRUE(v2 != nullptr);
const Vec3 grad2 = (*v2)(X);
EXPECT_DOUBLE_EQ((*v)(X),8.0);
EXPECT_NEAR(v->deriv(X,1),3.0,1e-6);
EXPECT_NEAR(grad1[0],3.0,1e-6);
EXPECT_NEAR(grad2[0],3.0,1e-6);
EXPECT_NEAR(v->deriv(X,2),4.0,1e-6);
EXPECT_NEAR(grad1[1],4.0,1e-6);
EXPECT_NEAR(grad2[1],4.0,1e-6);
EXPECT_NEAR(v->deriv(X,3),1.0,1e-6);
EXPECT_NEAR(grad1[2],1.0,1e-6);
EXPECT_NEAR(grad2[2],1.0,1e-6);
}
TEST(TestAnaSol, ParseADStress)
{
const char* input = "<anasol type=\"expression\" symmetric=\"true\" autodiff=\"true\">"
" <primary>"
" x^3+y^2+z | x^2+y^3+z"
" </primary>"
"</anasol>";
TiXmlDocument doc;
doc.Parse(input,nullptr,TIXML_ENCODING_UTF8);
const TiXmlElement* tag = doc.RootElement();
ASSERT_TRUE(tag != nullptr);
ASSERT_EQ(strcmp(tag->Value(),"anasol"),0);
AnaSol mySol(tag,false);
ASSERT_TRUE(mySol.getScalarSol() == nullptr);
ASSERT_TRUE(mySol.getVectorSol() != nullptr);
mySol.setupSecondarySolutions();
ASSERT_TRUE(mySol.getVectorSecSol() == nullptr);
const STensorFunc* v = mySol.getStressSol();
ASSERT_TRUE(v != nullptr);
Vec3 X(1.0,2.0,3.0);
const SymmTensor grad1 = (*v)(X);
EXPECT_DOUBLE_EQ(grad1(1,1), 3*X.x*X.x);
EXPECT_DOUBLE_EQ(grad1(1,2), 2*X.y);
EXPECT_DOUBLE_EQ(grad1(2,1), 2*X.y);
EXPECT_DOUBLE_EQ(grad1(2,2), 3*X.y*X.y);
}