From b996959c61de13c4023962cb6151ffca7247cc9a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 23 Oct 2013 20:09:55 +0200 Subject: [PATCH 1/5] Add new AutoDiffBlock::constant() overload without block sizes. This should simplify some uses of the autodiff code. The internals have been changed to allow for objects to have an empty vector of Jacobians, always treating that object as a constant. --- opm/autodiff/AutoDiffBlock.hpp | 57 +++++++++++++++++++++++++++++----- 1 file changed, 50 insertions(+), 7 deletions(-) diff --git a/opm/autodiff/AutoDiffBlock.hpp b/opm/autodiff/AutoDiffBlock.hpp index e9b6396c9..6e40e893a 100644 --- a/opm/autodiff/AutoDiffBlock.hpp +++ b/opm/autodiff/AutoDiffBlock.hpp @@ -98,6 +98,16 @@ namespace Opm /// Create an AutoDiffBlock representing a constant. /// \param[in] val values + static AutoDiffBlock constant(const V& val) + { + return AutoDiffBlock(val); + } + + /// Create an AutoDiffBlock representing a constant. + /// This variant allows to specify the block sizes used + /// for the Jacobian even though the Jacobian matrices + /// themselves will be zero. + /// \param[in] val values /// \param[in] blocksizes block pattern static AutoDiffBlock constant(const V& val, const std::vector& blocksizes) { @@ -165,14 +175,18 @@ namespace Opm /// Elementwise operator += AutoDiffBlock& operator+=(const AutoDiffBlock& rhs) { - assert (numBlocks() == rhs.numBlocks()); - assert (value().size() == rhs.value().size()); + if (jac_.empty()) { + jac_ = rhs.jac_; + } else if (!rhs.jac_.empty()) { + assert (numBlocks() == rhs.numBlocks()); + assert (value().size() == rhs.value().size()); - const int num_blocks = numBlocks(); - for (int block = 0; block < num_blocks; ++block) { - assert(jac_[block].rows() == rhs.jac_[block].rows()); - assert(jac_[block].cols() == rhs.jac_[block].cols()); - jac_[block] += rhs.jac_[block]; + const int num_blocks = numBlocks(); + for (int block = 0; block < num_blocks; ++block) { + assert(jac_[block].rows() == rhs.jac_[block].rows()); + assert(jac_[block].cols() == rhs.jac_[block].cols()); + jac_[block] += rhs.jac_[block]; + } } val_ += rhs.val_; @@ -183,6 +197,12 @@ namespace Opm /// Elementwise operator + AutoDiffBlock operator+(const AutoDiffBlock& rhs) const { + if (jac_.empty()) { + return val_ + rhs; + } + if (rhs.jac_.empty()) { + return *this - rhs.val_; + } std::vector jac = jac_; assert(numBlocks() == rhs.numBlocks()); int num_blocks = numBlocks(); @@ -197,6 +217,12 @@ namespace Opm /// Elementwise operator - AutoDiffBlock operator-(const AutoDiffBlock& rhs) const { + if (jac_.empty()) { + return val_ - rhs; + } + if (rhs.jac_.empty()) { + return *this - rhs.val_; + } std::vector jac = jac_; assert(numBlocks() == rhs.numBlocks()); int num_blocks = numBlocks(); @@ -211,6 +237,12 @@ namespace Opm /// Elementwise operator * AutoDiffBlock operator*(const AutoDiffBlock& rhs) const { + if (jac_.empty()) { + return val_ * rhs; + } + if (rhs.jac_.empty()) { + return *this * rhs.val_; + } int num_blocks = numBlocks(); std::vector jac(num_blocks); assert(numBlocks() == rhs.numBlocks()); @@ -228,6 +260,12 @@ namespace Opm /// Elementwise operator / AutoDiffBlock operator/(const AutoDiffBlock& rhs) const { + if (jac_.empty()) { + return val_ / rhs; + } + if (rhs.jac_.empty()) { + return *this / rhs.val_; + } int num_blocks = numBlocks(); std::vector jac(num_blocks); assert(numBlocks() == rhs.numBlocks()); @@ -292,6 +330,11 @@ namespace Opm } private: + AutoDiffBlock(const V& val) + : val_(val) + { + } + AutoDiffBlock(const V& val, const std::vector& jac) : val_(val), jac_(jac) From 46a17945a3176eb360c22fa647f02b1a81de8a67 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 23 Oct 2013 22:42:32 +0200 Subject: [PATCH 2/5] Fix minor whitespace issue. --- opm/autodiff/AutoDiffBlock.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/autodiff/AutoDiffBlock.hpp b/opm/autodiff/AutoDiffBlock.hpp index 6e40e893a..31175c49f 100644 --- a/opm/autodiff/AutoDiffBlock.hpp +++ b/opm/autodiff/AutoDiffBlock.hpp @@ -200,7 +200,7 @@ namespace Opm if (jac_.empty()) { return val_ + rhs; } - if (rhs.jac_.empty()) { + if (rhs.jac_.empty()) { return *this - rhs.val_; } std::vector jac = jac_; From 1eec8b16d6b79684dbbaa8b5c294a8827364d440 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 24 Oct 2013 00:14:32 +0200 Subject: [PATCH 3/5] Fix bug in operator+ introduced by previous commit. Also make documentation clearer. --- opm/autodiff/AutoDiffBlock.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/opm/autodiff/AutoDiffBlock.hpp b/opm/autodiff/AutoDiffBlock.hpp index 31175c49f..ad907dfd0 100644 --- a/opm/autodiff/AutoDiffBlock.hpp +++ b/opm/autodiff/AutoDiffBlock.hpp @@ -104,8 +104,8 @@ namespace Opm } /// Create an AutoDiffBlock representing a constant. - /// This variant allows to specify the block sizes used - /// for the Jacobian even though the Jacobian matrices + /// This variant requires specifying the block sizes used + /// for the Jacobians even though the Jacobian matrices /// themselves will be zero. /// \param[in] val values /// \param[in] blocksizes block pattern @@ -201,7 +201,7 @@ namespace Opm return val_ + rhs; } if (rhs.jac_.empty()) { - return *this - rhs.val_; + return *this + rhs.val_; } std::vector jac = jac_; assert(numBlocks() == rhs.numBlocks()); From 258d8e0e24b3d9c6b66f2572cf085bc3d74b3119 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 24 Oct 2013 13:41:55 +0200 Subject: [PATCH 4/5] Avoid infinite loop with two constant operands. --- opm/autodiff/AutoDiffBlock.hpp | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/opm/autodiff/AutoDiffBlock.hpp b/opm/autodiff/AutoDiffBlock.hpp index ad907dfd0..927927826 100644 --- a/opm/autodiff/AutoDiffBlock.hpp +++ b/opm/autodiff/AutoDiffBlock.hpp @@ -197,6 +197,9 @@ namespace Opm /// Elementwise operator + AutoDiffBlock operator+(const AutoDiffBlock& rhs) const { + if (jac_.empty() && rhs.jac_.empty()) { + return constant(val_ + rhs.val_); + } if (jac_.empty()) { return val_ + rhs; } @@ -217,6 +220,9 @@ namespace Opm /// Elementwise operator - AutoDiffBlock operator-(const AutoDiffBlock& rhs) const { + if (jac_.empty() && rhs.jac_.empty()) { + return constant(val_ - rhs.val_); + } if (jac_.empty()) { return val_ - rhs; } @@ -237,6 +243,9 @@ namespace Opm /// Elementwise operator * AutoDiffBlock operator*(const AutoDiffBlock& rhs) const { + if (jac_.empty() && rhs.jac_.empty()) { + return constant(val_ * rhs.val_); + } if (jac_.empty()) { return val_ * rhs; } @@ -260,6 +269,9 @@ namespace Opm /// Elementwise operator / AutoDiffBlock operator/(const AutoDiffBlock& rhs) const { + if (jac_.empty() && rhs.jac_.empty()) { + return constant(val_ / rhs.val_); + } if (jac_.empty()) { return val_ / rhs; } From 478c87dd749771b55d5481ef1a40b38eeef5cd25 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 25 Oct 2013 14:56:54 +0200 Subject: [PATCH 5/5] Reduce version to 0.9. Since there is a fair chance of interface adjustments after the release, the version number cautions the user a little bit more than calling it 1.0. --- dune.module | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dune.module b/dune.module index 96ffd60b0..9f59f898f 100644 --- a/dune.module +++ b/dune.module @@ -1,6 +1,6 @@ Module: opm-autodiff Description: Utilities for automatic differentiation and simulators based on AD -Version: 1.0 +Version: 0.9 Label: 2013.10 Maintainer: atgeirr@sintef.no Depends: opm-core