Commit 8c38539f authored by Oliver Sander's avatar Oliver Sander
Browse files

Implement a test for the subDiff method

parent 6dcb71c4
......@@ -19,6 +19,8 @@
#include <dune/istl/matrixindexset.hh>
#include <dune/solvers/common/interval.hh>
namespace Dune {
namespace TNNMG {
......@@ -247,6 +249,53 @@ void testDirectionalSubdifferential(const Functional& functional,
}
/** \brief Test the subDiff method
*
* This method tests the Functional::subDiff method by comparing its result with a
* Finite Difference approximation.
*
* \tparam Functional Type of the functional to be tested
* \param functional The functional to be tested
* \param testPoints Test subDiff at these points
*/
template <class Functional>
void testSubDiff(Functional& functional,
const std::vector<typename Functional::VectorType>& testPoints)
{
for (auto&& testPoint : testPoints)
{
// Step size. Best value: square root of the machine precision
const double eps = std::sqrt(std::numeric_limits<double>::epsilon());
for (size_t i=0; i<testPoint.size(); i++)
for (size_t j=0; j<testPoint[i].size(); j++)
{
Solvers::Interval<double> subDifferential;
// subDiff needs to be given the current point internally
functional.setVector(testPoint);
functional.subDiff(i, testPoint[i][j], subDifferential, j);
auto forwardPoint = testPoint;
forwardPoint[i][j] += eps;
auto backwardPoint = testPoint;
backwardPoint[i][j] -= eps;
auto forwardFDGradient = (functional(forwardPoint) - functional(testPoint)) / eps;
auto backwardFDGradient = (functional(testPoint) - functional(backwardPoint)) / eps;
if (std::abs(backwardFDGradient - subDifferential[0]) > 100*eps
or std::abs(forwardFDGradient - subDifferential[1]) > 100*eps)
{
std::cout << "Bug in subDiff for functional type " << Dune::className<Functional>() << ":" << std::endl;
std::cerr << "Subdifferential doesn't match FD approximation at coefficient (" << i << ", " << j << ")" << std::endl;
std::cerr << "SubDifferential: " << subDifferential
<< ", FD: [" << backwardFDGradient << ", " << forwardFDGradient << "]" << std::endl;
DUNE_THROW(MathError,"");
}
}
}
}
} // namespace TNNMG
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment