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 diff --git a/opm/autodiff/AutoDiffBlock.hpp b/opm/autodiff/AutoDiffBlock.hpp index e9b6396c9..927927826 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 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 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,15 @@ 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; + } + if (rhs.jac_.empty()) { + return *this + rhs.val_; + } std::vector jac = jac_; assert(numBlocks() == rhs.numBlocks()); int num_blocks = numBlocks(); @@ -197,6 +220,15 @@ 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; + } + if (rhs.jac_.empty()) { + return *this - rhs.val_; + } std::vector jac = jac_; assert(numBlocks() == rhs.numBlocks()); int num_blocks = numBlocks(); @@ -211,6 +243,15 @@ 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; + } + if (rhs.jac_.empty()) { + return *this * rhs.val_; + } int num_blocks = numBlocks(); std::vector jac(num_blocks); assert(numBlocks() == rhs.numBlocks()); @@ -228,6 +269,15 @@ 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; + } + if (rhs.jac_.empty()) { + return *this / rhs.val_; + } int num_blocks = numBlocks(); std::vector jac(num_blocks); assert(numBlocks() == rhs.numBlocks()); @@ -292,6 +342,11 @@ namespace Opm } private: + AutoDiffBlock(const V& val) + : val_(val) + { + } + AutoDiffBlock(const V& val, const std::vector& jac) : val_(val), jac_(jac)