diff --git a/AutoDiffBlock.hpp b/AutoDiffBlock.hpp index 4a8fe2706..fcf136ed6 100644 --- a/AutoDiffBlock.hpp +++ b/AutoDiffBlock.hpp @@ -156,6 +156,17 @@ namespace AutoDiff return jac_.size(); } + /// Sizes (number of columns) of Jacobian blocks. + std::vector blockPattern() const + { + const int nb = numBlocks(); + std::vector bp(nb); + for (int block = 0; block < nb; ++block) { + bp[block] = jac_[block].cols(); + } + return bp; + } + /// Function value const V& value() const { @@ -194,6 +205,7 @@ namespace AutoDiff return fw.print(os); } + /// Multiply with sparse matrix from the left. template ForwardBlock operator*(const typename ForwardBlock::M& lhs, @@ -210,6 +222,70 @@ namespace AutoDiff } + template + ForwardBlock operator*(const typename ForwardBlock::V& lhs, + const ForwardBlock& rhs) + { + return ForwardBlock::constant(lhs, rhs.blockPattern()) * rhs; + } + + + template + ForwardBlock operator*(const ForwardBlock& lhs, + const typename ForwardBlock::V& rhs) + { + return rhs * lhs; // Commutative operation. + } + + + template + ForwardBlock operator+(const typename ForwardBlock::V& lhs, + const ForwardBlock& rhs) + { + return ForwardBlock::constant(lhs, rhs.blockPattern()) + rhs; + } + + + template + ForwardBlock operator+(const ForwardBlock& lhs, + const typename ForwardBlock::V& rhs) + { + return rhs + lhs; // Commutative operation. + } + + + template + ForwardBlock operator-(const typename ForwardBlock::V& lhs, + const ForwardBlock& rhs) + { + return ForwardBlock::constant(lhs, rhs.blockPattern()) - rhs; + } + + + template + ForwardBlock operator-(const ForwardBlock& lhs, + const typename ForwardBlock::V& rhs) + { + return lhs - ForwardBlock::constant(rhs, lhs.blockPattern()); + } + + + template + ForwardBlock operator/(const typename ForwardBlock::V& lhs, + const ForwardBlock& rhs) + { + return ForwardBlock::constant(lhs, rhs.blockPattern()) / rhs; + } + + + template + ForwardBlock operator/(const ForwardBlock& lhs, + const typename ForwardBlock::V& rhs) + { + return lhs / ForwardBlock::constant(rhs, lhs.blockPattern()); + } + + } // namespace Autodiff diff --git a/AutoDiffBlockArma.hpp b/AutoDiffBlockArma.hpp new file mode 100644 index 000000000..91fda60d9 --- /dev/null +++ b/AutoDiffBlockArma.hpp @@ -0,0 +1,216 @@ +/* + Copyright 2013 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_AUTODIFFBLOCKARMA_HEADER_INCLUDED +#define OPM_AUTODIFFBLOCKARMA_HEADER_INCLUDED + + +// #include "AutoDiff.hpp" +// #include +// #include +#include +#include +#include + +namespace AutoDiff +{ + + template + class ForwardBlock + { + public: + /// Underlying types for scalar vectors and jacobians. + typedef arma::Col V; + typedef arma::SpMat M; + + /// Named constructor pattern used here. + static ForwardBlock constant(const V& val, const std::vector& blocksizes) + { + std::vector jac; + const int num_elem = val.size(); + const int num_blocks = blocksizes.size(); + // For constants, all jacobian blocks are zero. + jac.resize(num_blocks); + for (int i = 0; i < num_blocks; ++i) { + jac[i] = M(num_elem, blocksizes[i]); + } + return ForwardBlock(val, jac); + } + + static ForwardBlock variable(const int index, const V& val, const std::vector& blocksizes) + { + std::vector jac; + const int num_elem = val.size(); + const int num_blocks = blocksizes.size(); + // First, set all jacobian blocks to zero... + jac.resize(num_blocks); + for (int i = 0; i < num_blocks; ++i) { + jac[i] = M(num_elem, blocksizes[i]); + } + // ... then set the one corrresponding to this variable to identity. + assert(blocksizes[index] == num_elem); + jac[index].eye(num_elem, num_elem); + return ForwardBlock(val, jac); + } + + static ForwardBlock function(const V& val, const std::vector& jac) + { + return ForwardBlock(val, jac); + } + + /// Operator + + ForwardBlock operator+(const ForwardBlock& rhs) + { + std::vector jac = jac_; + assert(numBlocks() == rhs.numBlocks()); + int num_blocks = numBlocks(); + for (int block = 0; block < num_blocks; ++block) { + assert(jac[block].n_rows == rhs.jac_[block].n_rows); + assert(jac[block].n_cols == rhs.jac_[block].n_cols); + jac[block] += rhs.jac_[block]; + } + return function(val_ + rhs.val_, jac); + } + + /// Operator - + ForwardBlock operator-(const ForwardBlock& rhs) + { + std::vector jac = jac_; + assert(numBlocks() == rhs.numBlocks()); + int num_blocks = numBlocks(); + for (int block = 0; block < num_blocks; ++block) { + assert(jac[block].n_rows == rhs.jac_[block].n_rows); + assert(jac[block].n_cols == rhs.jac_[block].n_cols); + jac[block] -= rhs.jac_[block]; + } + return function(val_ - rhs.val_, jac); + } + + /// Operator * + ForwardBlock operator*(const ForwardBlock& rhs) + { + int num_blocks = numBlocks(); + std::vector jac(num_blocks); + assert(numBlocks() == rhs.numBlocks()); + typedef Eigen::DiagonalMatrix D; + D D1 = val_.matrix().asDiagonal(); + D D2 = rhs.val_.matrix().asDiagonal(); + for (int block = 0; block < num_blocks; ++block) { + assert(jac_[block].n_rows == rhs.jac_[block].n_rows); + assert(jac_[block].n_cols == rhs.jac_[block].n_cols); + jac[block] = D2*jac_[block] + D1*rhs.jac_[block]; + } + return function(val_ * rhs.val_, jac); + } + + /// Operator / + ForwardBlock operator/(const ForwardBlock& rhs) + { + int num_blocks = numBlocks(); + std::vector jac(num_blocks); + assert(numBlocks() == rhs.numBlocks()); + typedef Eigen::DiagonalMatrix D; + D D1 = val_.matrix().asDiagonal(); + D D2 = rhs.val_.matrix().asDiagonal(); + D D3 = std::pow(rhs.val_, -2).matrix().asDiagonal(); + for (int block = 0; block < num_blocks; ++block) { + assert(jac_[block].n_rows == rhs.jac_[block].n_rows); + assert(jac_[block].n_cols == rhs.jac_[block].n_cols); + jac[block] = D3 * (D2*jac_[block] - D1*rhs.jac_[block]); + } + return function(val_ / rhs.val_, jac); + } + /// I/O. + template + Ostream& + print(Ostream& os) const + { + int num_blocks = jac_.size(); + os << "Value =\n" << val_ << "\n\nJacobian =\n"; + for (int i = 0; i < num_blocks; ++i) { + os << "Sub Jacobian #" << i << '\n' << jac_[i] << "\n"; + } + return os; + } + + /// Number of variables or Jacobian blocks. + int numBlocks() const + { + return jac_.size(); + } + + /// Function value + const V& value() const + { + return val_; + } + + /// Function derivatives + const std::vector& derivative() const + { + return jac_; + } + + private: + ForwardBlock(const V& val, + const std::vector& jac) + : val_(val), jac_(jac) + { +#ifndef NDEBUG + const int num_elem = val_.size(); + const int num_blocks = jac_.size(); + for (int block = 0; block < num_blocks; ++block) { + assert(num_elem == jac_[block].n_rows); + } +#endif + } + + V val_; + std::vector jac_; + }; + + + template + Ostream& + operator<<(Ostream& os, const ForwardBlock& fw) + { + return fw.print(os); + } + + /// Multiply with sparse matrix from the left. + template + ForwardBlock operator*(const typename ForwardBlock::M& lhs, + const ForwardBlock& rhs) + { + int num_blocks = rhs.numBlocks(); + std::vector::M> jac(num_blocks); + assert(lhs.n_cols == rhs.value().n_rows); + for (int block = 0; block < num_blocks; ++block) { + jac[block] = lhs*rhs.derivative()[block]; + } + typename ForwardBlock::V val = lhs*rhs.value().matrix(); + return ForwardBlock::function(val, jac); + } + + +} // namespace Autodiff + + + +#endif // OPM_AUTODIFFBLOCKARMA_HEADER_INCLUDED diff --git a/sim_simple.cpp b/sim_simple.cpp index aced416bf..59b1b3608 100644 --- a/sim_simple.cpp +++ b/sim_simple.cpp @@ -55,8 +55,8 @@ struct HelperOps typedef AutoDiff::ForwardBlock::V V; /// A list of internal faces. - typedef Eigen::Array iFaces; - iFaces internal_faces; + typedef Eigen::Array IFaces; + IFaces internal_faces; /// Extract for each face the difference of its adjacent cells'values. M ngrad; @@ -172,13 +172,13 @@ namespace { select(const typename ADB::V& press, const std::vector& xc ) const { - typedef HelperOps::iFaces::Index ifIndex; - const ifIndex nif = h_.internal_faces.size(); + typedef HelperOps::IFaces::Index IFIndex; + const IFIndex nif = h_.internal_faces.size(); // Define selector structure. typedef typename Eigen::Triplet Triplet; std::vector s; s.reserve(nif); - for (ifIndex i = 0; i < nif; ++i) { + for (IFIndex i = 0; i < nif; ++i) { const int f = h_.internal_faces[i]; const int c1 = g_.face_cells[2*f + 0]; const int c2 = g_.face_cells[2*f + 1]; @@ -276,27 +276,27 @@ int main() Opm::time::StopWatch clock; clock.start(); - Opm::GridManager gm(3,3);//(50, 50, 10); + const Opm::GridManager gm(3,3);//(50, 50, 10); const UnstructuredGrid& grid = *gm.c_grid(); using namespace Opm::unit; using namespace Opm::prefix; - // Opm::IncompPropertiesBasic props(2, Opm::SaturationPropsBasic::Linear, - // { 1000.0, 800.0 }, - // { 1.0*centi*Poise, 5.0*centi*Poise }, - // 0.2, 100*milli*darcy, - // grid.dimensions, grid.number_of_cells); - // Opm::IncompPropertiesBasic props(2, Opm::SaturationPropsBasic::Linear, - // { 1000.0, 1000.0 }, - // { 1.0, 1.0 }, - // 1.0, 1.0, - // grid.dimensions, grid.number_of_cells); - Opm::IncompPropertiesBasic props(2, Opm::SaturationPropsBasic::Linear, - { 1000.0, 1000.0 }, - { 1.0, 30.0 }, - 1.0, 1.0, - grid.dimensions, grid.number_of_cells); + // const Opm::IncompPropertiesBasic props(2, Opm::SaturationPropsBasic::Linear, + // { 1000.0, 800.0 }, + // { 1.0*centi*Poise, 5.0*centi*Poise }, + // 0.2, 100*milli*darcy, + // grid.dimensions, grid.number_of_cells); + // const Opm::IncompPropertiesBasic props(2, Opm::SaturationPropsBasic::Linear, + // { 1000.0, 1000.0 }, + // { 1.0, 1.0 }, + // 1.0, 1.0, + // grid.dimensions, grid.number_of_cells); + const Opm::IncompPropertiesBasic props(2, Opm::SaturationPropsBasic::Linear, + { 1000.0, 1000.0 }, + { 1.0, 30.0 }, + 1.0, 1.0, + grid.dimensions, grid.number_of_cells); std::vector htrans(grid.cell_facepos[grid.number_of_cells]); - tpfa_htrans_compute((UnstructuredGrid*)&grid, props.permeability(), htrans.data()); + tpfa_htrans_compute(const_cast(&grid), props.permeability(), htrans.data()); // std::vector trans(grid.number_of_faces); V trans_all(grid.number_of_faces); tpfa_trans_compute((UnstructuredGrid*)&grid, htrans.data(), trans_all.data()); @@ -308,7 +308,7 @@ int main() std::cerr << "Opm core " << clock.secsSinceLast() << std::endl; // Define neighbourhood-derived operator matrices. - HelperOps ops(grid); + const HelperOps ops(grid); const int num_internal = ops.internal_faces.size(); V transi(num_internal); for (int fi = 0; fi < num_internal; ++fi) { @@ -334,15 +334,15 @@ int main() // totmob - explicit as well TwoCol kr(nc, 2); props.relperm(nc, s0.data(), allcells.data(), kr.data(), 0); - V krw = kr.leftCols<1>(); - V kro = kr.rightCols<1>(); + const V krw = kr.leftCols<1>(); + const V kro = kr.rightCols<1>(); const double* mu = props.viscosity(); - V totmob = krw/mu[0] + kro/mu[1]; - V totmobf = (ops.caver*totmob.matrix()).array(); + const V totmob = krw/mu[0] + kro/mu[1]; + const V totmobf = (ops.caver*totmob.matrix()).array(); // Mobility-weighted transmissibilities per internal face. // Still explicit, and no upwinding! - V mobtransf = totmobf*transi; + const V mobtransf = totmobf*transi; std::cerr << "Property arrays " << clock.secsSinceLast() << std::endl; @@ -351,18 +351,14 @@ int main() p0.fill(200*Opm::unit::barsa); // First actual AD usage: defining pressure variable. - std::vector block_pattern = { nc }; - // Could actually write { nc } instead of block_pattern below, + const std::vector bpat = { nc }; + // Could actually write { nc } instead of bpat below, // but we prefer a named variable since we will repeat it. - ADB p = ADB::variable(0, p0, block_pattern); - ADB ngradp = ops.ngrad*p; + const ADB p = ADB::variable(0, p0, bpat); + const ADB ngradp = ops.ngrad*p; // We want flux = totmob*trans*(p_i - p_j) for the ij-face. - // We only need to multiply mobtransf and pdiff_face, - // but currently multiplication with constants is not in, - // so we define an AD constant to multiply with. - ADB mobtransf_ad = ADB::constant(mobtransf, block_pattern); - ADB flux = mobtransf_ad*ngradp; - ADB residual = ops.div*flux - ADB::constant(q, block_pattern); + const ADB flux = mobtransf*ngradp; + const ADB residual = ops.div*flux - q; std::cerr << "Construct AD residual " << clock.secsSinceLast() << std::endl; // It's the residual we want to be zero. We know it's linear in p, @@ -374,20 +370,21 @@ int main() // residual.derived()[0]. Eigen::UmfPackLU solver; - M matr = residual.derivative()[0]; - matr.coeffRef(0,0) *= 2.0; - matr.makeCompressed(); - solver.compute(matr); + M pmatr = residual.derivative()[0]; + pmatr.coeffRef(0,0) *= 2.0; + pmatr.makeCompressed(); + solver.compute(pmatr); if (solver.info() != Eigen::Success) { std::cerr << "Pressure/flow Jacobian decomposition error\n"; return EXIT_FAILURE; } - Eigen::VectorXd x = solver.solve(residual.value().matrix()); + // const Eigen::VectorXd dp = solver.solve(residual.value().matrix()); + const V dp = solver.solve(residual.value().matrix()).array(); if (solver.info() != Eigen::Success) { std::cerr << "Pressure/flow solve failure\n"; return EXIT_FAILURE; } - V p1 = p0 - x.array(); + const V p1 = p0 - dp; std::cerr << "Solve " << clock.secsSinceLast() << std::endl; // std::cout << p1 << std::endl; @@ -401,69 +398,52 @@ int main() // and f_w is (for now) based on averaged mobilities, not upwind. double res_norm = 1e100; - V s1 = /*s0.leftCols<1>()*/0.5*V::Ones(nc,1); // Initial guess. - UpwindSelector upws(grid, ops); - const ADB nkdp = (ADB::constant(transi , block_pattern) * - ADB::constant(ops.ngrad * p1.matrix(), block_pattern)); - const ADB s00 = ADB::constant(s0.leftCols<1>(), block_pattern); - const std::vector pmobc0 = phaseMobility(props, allcells, s00.value()); - const std::vector pmobf0 = upws.select(p1, pmobc0); - const std::vector null = { ADB::M(transi.size(), nc) }; - const ADB dflux = (ADB::function((pmobf0[0] + pmobf0[1]).value(), null) * - ADB::function(nkdp.value() , null)); + const V sw0 = s0.leftCols<1>(); + // V sw1 = sw0; + V sw1 = 0.5*V::Ones(nc,1); + const UpwindSelector upws(grid, ops); + const V nkdp = transi * (ops.ngrad * p1.matrix()).array(); + const V dflux = totmobf * nkdp; + const V pv = Eigen::Map(props.porosity(), nc, 1) + * Eigen::Map(grid.cell_volumes, nc, 1); + const double dt = 0.0005; + const V dtpv = dt/pv; + const V qneg = q.min(V::Zero(nc,1)); + const V qpos = q.max(V::Zero(nc,1)); std::cout.setf(std::ios::scientific); std::cout.precision(16); int it = 0; do { - const std::vector& bp = block_pattern; - ADB s = ADB::variable(0, s1, bp); - const double dt = 0.0005; - V pv = Eigen::Map(props.porosity(), nc, 1) - * Eigen::Map(grid.cell_volumes, nc, 1); - V dtpv = dt/pv; - // std::cout << dtpv; - std::vector pmobc = phaseMobility(props, allcells, s.value()); - std::vector pmobf = upws.select(p1, pmobc); - - ADB fw_cell = fluxFunc(pmobc); + const ADB sw = ADB::variable(0, sw1, bpat); + const std::vector pmobc = phaseMobility(props, allcells, sw.value()); + const std::vector pmobf = upws.select(p1, pmobc); + const ADB fw_cell = fluxFunc(pmobc); const ADB fw_face = fluxFunc(pmobf); - ADB flux1 = fw_face * dflux; - // std::cout << flux1; - V qneg = dtpv*q; - V qpos = dtpv*q; - // Cheating a bit... - qneg[0] = 0.0; - qpos[nc-1] = 0.0; - ADB qtr_ad = ADB::constant(qpos, bp) + fw_cell*ADB::constant(qneg, bp); - ADB transport_residual = s - ADB::constant(s0.leftCols<1>(), bp) - + ADB::constant(dtpv, bp)*(ops.div*flux1) - - qtr_ad; + const ADB flux1 = fw_face * dflux; + const ADB qtr_ad = qpos + fw_cell*qneg; + const ADB transport_residual = sw - sw0 + dtpv*(ops.div*flux1 - qtr_ad); res_norm = transport_residual.value().matrix().norm(); std::cout << "res_norm[" << it << "] = " << res_norm << std::endl; - matr = transport_residual.derivative()[0]; - matr.makeCompressed(); - // std::cout << transport_residual; - solver.compute(matr); + M smatr = transport_residual.derivative()[0]; + smatr.makeCompressed(); + solver.compute(smatr); if (solver.info() != Eigen::Success) { std::cerr << "Transport Jacobian decomposition error\n"; return EXIT_FAILURE; } - x = solver.solve(transport_residual.value().matrix()); + const V ds = solver.solve(transport_residual.value().matrix()).array(); if (solver.info() != Eigen::Success) { std::cerr << "Transport solve failure\n"; return EXIT_FAILURE; } - // std::cout << x << std::endl; - s1 = s.value() - x.array(); + sw1 = sw.value() - ds; std::cerr << "Solve for s[" << it << "]: " << clock.secsSinceLast() << '\n'; - for (int c = 0; c < nc; ++c) { - s1[c] = std::min(1.0, std::max(0.0, s1[c])); - } + sw1 = sw1.min(V::Ones(nc,1)).max(V::Zero(nc,1)); it += 1; } while (res_norm > 1e-7);