diff --git a/sim_simple.cpp b/sim_simple.cpp index fa6cd299f..01dc14063 100644 --- a/sim_simple.cpp +++ b/sim_simple.cpp @@ -44,6 +44,71 @@ */ +/// Contains vectors and sparse matrices that represent subsets or +/// operations on (AD or regular) vectors of data. +struct HelperOps +{ + typedef AutoDiff::ForwardBlock::M M; + typedef AutoDiff::ForwardBlock::V V; + + /// A list of internal faces. + Eigen::Array internal_faces; + + /// Extract for each face the difference of its adjacent cells'values. + M ngrad; + /// Extract for each face the average of its adjacent cells' values. + M caver; + /// Extract for each cell the sum of its adjacent faces' (signed) values. + M div; + + /// Constructs all helper vectors and matrices. + HelperOps(const UnstructuredGrid& grid) + { + const int nc = grid.number_of_cells; + const int nf = grid.number_of_faces; + // Define some neighbourhood-derived helper arrays. + typedef Eigen::Array OneColInt; + typedef Eigen::Array OneColBool; + typedef Eigen::Array TwoColInt; + typedef Eigen::Array TwoColBool; + TwoColInt nb = Eigen::Map(grid.face_cells, nf, 2); + // std::cout << "nb = \n" << nb << std::endl; + TwoColBool nbib = nb >= 0; + OneColBool ifaces = nbib.rowwise().all(); + const int num_internal = ifaces.cast().sum(); + // std::cout << num_internal << " internal faces." << std::endl; + TwoColInt nbi(num_internal, 2); + internal_faces.resize(num_internal); + int fi = 0; + for (int f = 0; f < nf; ++f) { + if (ifaces[f]) { + internal_faces[fi] = f; + nbi.row(fi) = nb.row(f); + ++fi; + } + } + // std::cout << "nbi = \n" << nbi << std::endl; + // Create matrices. + ngrad.resize(num_internal, nc); + caver.resize(num_internal, nc); + typedef Eigen::Triplet Tri; + std::vector ngrad_tri; + std::vector caver_tri; + ngrad_tri.reserve(2*num_internal); + caver_tri.reserve(2*num_internal); + for (int i = 0; i < num_internal; ++i) { + ngrad_tri.emplace_back(i, nbi(i,0), 1.0); + ngrad_tri.emplace_back(i, nbi(i,1), -1.0); + caver_tri.emplace_back(i, nbi(i,0), 0.5); + caver_tri.emplace_back(i, nbi(i,1), 0.5); + } + ngrad.setFromTriplets(ngrad_tri.begin(), ngrad_tri.end()); + caver.setFromTriplets(caver_tri.begin(), caver_tri.end()); + div = ngrad.transpose(); + } +}; + + int main() { @@ -68,55 +133,19 @@ int main() V trans_all(grid.number_of_faces); tpfa_trans_compute((UnstructuredGrid*)&grid, htrans.data(), trans_all.data()); const int nc = grid.number_of_cells; - const int nf = grid.number_of_faces; std::vector allcells(nc); for (int i = 0; i < nc; ++i) { allcells[i] = i; } std::cerr << "Opm core " << clock.secsSinceLast() << std::endl; - // Define neighbourhood-derived matrices. - typedef Eigen::Array OneColInt; - typedef Eigen::Array OneColBool; - typedef Eigen::Array TwoColInt; - typedef Eigen::Array TwoColBool; - TwoColInt nb = Eigen::Map(grid.face_cells, nf, 2); - // std::cout << "nb = \n" << nb << std::endl; - TwoColBool nbib = nb >= 0; - OneColBool ifaces = nbib.rowwise().all(); - const int num_internal = ifaces.cast().sum(); - // std::cout << num_internal << " internal faces." << std::endl; - TwoColInt nbi(num_internal, 2); + // Define neighbourhood-derived operator matrices. + HelperOps ops(grid); + const int num_internal = ops.internal_faces.size(); V transi(num_internal); - int fi = 0; - for (int f = 0; f < nf; ++f) { - if (ifaces[f]) { - transi[fi] = trans_all[f]; - nbi.row(fi) = nb.row(f); - ++fi; - } + for (int fi = 0; fi < num_internal; ++fi) { + transi[fi] = trans_all[ops.internal_faces[fi]]; } - // std::cout << "nbi = \n" << nbi << std::endl; - // Create matrices: - // cdiff - a matrix for computing cell-cell differences per face. - // caver - a matrix for computing cell-cell averages per face. - // div - a matrix for computing divergence at a cell from face-given fluxes. - M cdiff(num_internal, nc); - M caver(num_internal, nc); - typedef Eigen::Triplet Tri; - std::vector cdiff_tri; - std::vector caver_tri; - cdiff_tri.reserve(2*num_internal); - caver_tri.reserve(2*num_internal); - for (int i = 0; i < num_internal; ++i) { - cdiff_tri.emplace_back(i, nbi(i,0), 1.0); - cdiff_tri.emplace_back(i, nbi(i,1), -1.0); - caver_tri.emplace_back(i, nbi(i,0), 0.5); - caver_tri.emplace_back(i, nbi(i,1), 0.5); - } - cdiff.setFromTriplets(cdiff_tri.begin(), cdiff_tri.end()); - caver.setFromTriplets(caver_tri.begin(), caver_tri.end()); - M div = cdiff.transpose(); std::cerr << "Topology matrices " << clock.secsSinceLast() << std::endl; typedef AutoDiff::ForwardBlock ADB; @@ -141,7 +170,7 @@ int main() V kro = kr.rightCols<1>(); const double* mu = props.viscosity(); V totmob = krw/mu[0] + kro/mu[1]; - V totmobf = (caver*totmob.matrix()).array(); + V totmobf = (ops.caver*totmob.matrix()).array(); // Mobility-weighted transmissibilities per internal face. // Still explicit, and no upwinding! @@ -149,20 +178,23 @@ int main() std::cerr << "Property arrays " << clock.secsSinceLast() << std::endl; + // Initial pressure. + V p0(nc,1); + p0.fill(200*Opm::unit::barsa); - // First actual AD usage: defining pressure. + // First actual AD usage: defining pressure variable. std::vector block_pattern = { nc }; // Could actually write { nc } instead of block_pattern below, // but we prefer a named variable since we will repeat it. - ADB p = ADB::variable(0, V::Zero(nc, 1), block_pattern); - ADB pdiff_face = cdiff*p; + ADB p = ADB::variable(0, p0, block_pattern); + 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*pdiff_face; - ADB residual = div*flux - ADB::constant(q, block_pattern); + ADB flux = mobtransf_ad*ngradp; + ADB residual = ops.div*flux - ADB::constant(q, block_pattern); std::cerr << "Construct AD residual " << clock.secsSinceLast() << std::endl; // std::cout << div << pdiff_face; @@ -171,7 +203,12 @@ int main() // std::cout << residual << std::endl; // It's the residual we want to be zero. We know it's linear in p, - // so we just need a single linear solve. + // so we just need a single linear solve. Since we have formulated + // ourselves with a residual and jacobian we do this with a single + // Newton step (hopefully easy to extend later): + // p = p0 - J(p0) \ R(p0) + // Where R(p0) and J(p0) are contained in residual.value() and + // residual.derived()[0]. Eigen::UmfPackLU solver; M matr = residual.derivative()[0]; @@ -187,6 +224,7 @@ int main() // std::cerr << "Solve failure!\n"; // return 1; // } + V p_new = p0 - x.array(); std::cerr << "Solve " << clock.secsSinceLast() << std::endl; - std::cout << x << std::endl; + std::cout << p_new << std::endl; }