diff --git a/sim_simple.cpp b/sim_simple.cpp index c889040db..5833e41ad 100644 --- a/sim_simple.cpp +++ b/sim_simple.cpp @@ -110,6 +110,43 @@ struct HelperOps } }; +/// Returns fw(sw). +template +ADB fluxFunc(const Opm::IncompPropertiesInterface& props, + const std::vector& cells, + const typename ADB::V& sw) +{ + typedef Eigen::Array TwoCol; + typedef Eigen::Array FourCol; + typedef typename ADB::V V; + typedef typename ADB::M M; + const int nc = props.numCells(); + TwoCol s(nc, 2); + s.leftCols<1>() = sw; + s.rightCols<1>() = 1.0 - s.leftCols<1>(); + TwoCol kr(nc, 2); + FourCol dkr(nc, 4); + props.relperm(nc, s.data(), cells.data(), kr.data(), dkr.data()); + V krw = kr.leftCols<1>(); + V kro = kr.rightCols<1>(); + V dkrw = dkr.leftCols<1>(); // Left column is top-left of dkr/ds 2x2 matrix. + V dkro = -dkr.rightCols<1>(); // Right column is bottom-right of dkr/ds 2x2 matrix. + M krwjac(nc,nc); + M krojac(nc,nc); + auto sizes = Eigen::ArrayXi::Ones(nc); + krwjac.reserve(sizes); + krojac.reserve(sizes); + for (int c = 0; c < nc; ++c) { + krwjac.insert(c,c) = dkrw(c); + krojac.insert(c,c) = dkro(c); + } + const double* mu = props.viscosity(); + ADB mw_ad = ADB::function(krw/mu[0], { krwjac/mu[0] }); + ADB mo_ad = ADB::function(kro/mu[1], { krojac/mu[1] }); + ADB fw = mw_ad / (mw_ad + mo_ad); + // std::cout << mw_ad << mo_ad << (mw_ad + mo_ad) << fw; + return fw; +} int main() @@ -120,14 +157,24 @@ int main() Opm::time::StopWatch clock; clock.start(); - Opm::GridManager gm(50, 50, 10); + 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::Quadratic, - { 1000.0, 800.0 }, - { 1.0*centi*Poise, 5.0*centi*Poise }, - 0.2, 100*milli*darcy, + // 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); std::vector htrans(grid.cell_facepos[grid.number_of_cells]); tpfa_htrans_compute((UnstructuredGrid*)&grid, props.permeability(), htrans.data()); @@ -159,15 +206,15 @@ int main() q[0] = 1.0; q[nc-1] = -1.0; - // s - this is explicit now + // s0 - this is explicit now typedef Eigen::Array TwoCol; - TwoCol s(nc, 2); - s.leftCols<1>().setZero(); - s.rightCols<1>().setOnes(); + TwoCol s0(nc, 2); + s0.leftCols<1>().setZero(); + s0.rightCols<1>().setOnes(); // totmob - explicit as well TwoCol kr(nc, 2); - props.relperm(nc, s.data(), allcells.data(), kr.data(), 0); + props.relperm(nc, s0.data(), allcells.data(), kr.data(), 0); V krw = kr.leftCols<1>(); V kro = kr.rightCols<1>(); const double* mu = props.viscosity(); @@ -199,11 +246,6 @@ int main() ADB residual = ops.div*flux - ADB::constant(q, block_pattern); std::cerr << "Construct AD residual " << clock.secsSinceLast() << std::endl; - // 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. Since we have formulated // ourselves with a residual and jacobian we do this with a single @@ -226,7 +268,68 @@ int main() // std::cerr << "Solve failure!\n"; // return 1; // } - V p_new = p0 - x.array(); + V p1 = p0 - x.array(); std::cerr << "Solve " << clock.secsSinceLast() << std::endl; - std::cout << p_new << std::endl; + // std::cout << p1 << std::endl; + + // ------ Transport solve ------ + + // Now we'll try to do a transport step as well. + // Residual formula is + // R_w = s_w - s_w^0 + dt/pv * (div v_w) + // where + // v_w = f_w v + // 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. + 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; + V ngradp1 = ops.ngrad*p1.matrix(); + // std::cout << ngradp1 << std::endl; + ADB fw_cell = fluxFunc(props, allcells, s.value()); + // std::cout << fw_cell; + ADB fw_face = ops.caver*fw_cell; + // std::cout << fw_face; + ADB flux1 = fw_face*ADB::constant(ngradp1, bp); + // 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; + res_norm = transport_residual.value().matrix().norm(); + std::cout << "res_norm = " << res_norm << std::endl; + + matr = transport_residual.derivative()[0]; + matr.makeCompressed(); + // std::cout << transport_residual; + solver.compute(matr); + // if (solver.info() != Eigen::Succeeded) { + // std::cerr << "Decomposition error!\n"; + // return 1; + // } + x = solver.solve(transport_residual.value().matrix()); + // if (solver.info() != Eigen::Succeeded) { + // std::cerr << "Solve failure!\n"; + // return 1; + // } + // std::cout << x << std::endl; + s1 = s.value() - x.array(); + std::cerr << "Solve for s " << clock.secsSinceLast() << std::endl; + for (int c = 0; c < nc; ++c) { + s1[c] = std::min(1.0, std::max(0.0, s1[c])); + } + std::cout << "s1 = \n" << s1 << std::endl; + } while (res_norm > 1e-7); }