From 169557111d756cf690f3351b0fcd58c57b36b7b3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 2 May 2013 08:51:36 +0200 Subject: [PATCH] Simple simulator-like test program added. So far it only computes a single pressure solution. --- sim_simple.cpp | 177 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 177 insertions(+) create mode 100644 sim_simple.cpp diff --git a/sim_simple.cpp b/sim_simple.cpp new file mode 100644 index 000000000..14c39ea52 --- /dev/null +++ b/sim_simple.cpp @@ -0,0 +1,177 @@ +/* + 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 . +*/ + +#include "AutoDiffBlock.hpp" +#include +#include +#include +#include +#include +#include + +/* + Equations for incompressible two-phase flow. + + Using s and p as variables: + + PV (s_i - s0_i) / dt + sum_{j \in U(i)} f(s_j) v_{ij} + sum_{j in D(i) f(s_i) v_{ij} = qw_i + + where + + v_{ij} = totmob_ij T_ij (p_i - p_j) + + + Pressure equation: + + sum_{j \in N(i)} totmob_ij T_ij (p_i - p_j) = q_i + +*/ + + +int main() +{ + typedef AutoDiff::ForwardBlock ADB; + typedef ADB::V V; + typedef ADB::M M; + + Opm::GridManager gm(100,100); + const UnstructuredGrid& grid = *gm.c_grid(); + using namespace Opm::unit; + using namespace Opm::prefix; + Opm::IncompPropertiesBasic props(2, Opm::SaturationPropsBasic::Quadratic, + { 1000.0, 800.0 }, + { 1.0*centi*Poise, 5.0*centi*Poise }, + 0.2, 100*milli*darcy, + 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()); + // std::vector trans(grid.number_of_faces); + 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; + } + // 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); + 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; + } + } + // 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); + for (int i = 0; i < num_internal; ++i) { + cdiff.insert(i, nbi(i,0)) = 1.0; + cdiff.insert(i, nbi(i,1)) = -1.0; + caver.insert(i, nbi(i,0)) = 0.5; + caver.insert(i, nbi(i,1)) = 0.5; + } + M div = cdiff.transpose(); + + typedef AutoDiff::ForwardBlock ADB; + typedef ADB::V V; + + // q + V q(nc); + q.setZero(); + q[0] = 1.0; + q[nc-1] = -1.0; + + // s - this is explicit now + typedef Eigen::Array TwoCol; + TwoCol s(nc, 2); + s.leftCols<1>().setZero(); + s.rightCols<1>().setOnes(); + + V sw = s.leftCols<1>(); + V so = s.rightCols<1>(); + + // totmob - explicit as well + TwoCol kr(nc, 2); + props.relperm(nc, s.data(), allcells.data(), kr.data(), 0); + V krw = kr.leftCols<1>(); + V kro = kr.rightCols<1>(); + const double* mu = props.viscosity(); + V totmob = krw/mu[0] + kro/mu[1]; + V totmobf = (caver*totmob.matrix()).array(); + + // Mobility-weighted transmissibilities per internal face. + // Still explicit, and no upwinding! + V mobtransf = totmobf*transi; + + // First actual AD usage: defining pressure. + 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; + // 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); + + // std::cout << div << pdiff_face; + // std::cout << div*pdiff_face; + // std::cout << q << std::endl; + // 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. + + Eigen::UmfPackLU solver; + M matr = residual.derivative()[0]; + matr.coeffRef(0,0) *= 2.0; + matr.makeCompressed(); + solver.compute(residual.derivative()[0]); + // if (solver.info() != Eigen::Succeeded) { + // std::cerr << "Decomposition error!\n"; + // return 1; + // } + Eigen::VectorXd x = solver.solve(residual.value().matrix()); + // if (solver.info() != Eigen::Succeeded) { + // std::cerr << "Solve failure!\n"; + // return 1; + // } + std::cout << x << std::endl; +}