First partitioning into source/examples/tests

Reorganisations are likely.

examples/find_zero.cpp

@@ -0,0 +1,102 @@
+/*===========================================================================
+//
+// File: find_zero.cpp
+//
+// Created: 2013-04-29 11:58:29+0200
+//
+// Authors: Knut-Andreas Lie      <Knut-Andreas.Lie@sintef.no>
+//          Halvor M. Nilsen      <HalvorMoll.Nilsen@sintef.no>
+//          Atgeirr F. Rasmussen  <atgeirr@sintef.no>
+//          Xavier Raynaud        <Xavier.Raynaud@sintef.no>
+//          Bård Skaflestad       <Bard.Skaflestad@sintef.no>
+//
+//==========================================================================*/
+
+
+/*
+  Copyright 2013 SINTEF ICT, Applied Mathematics.
+  Copyright 2013 Statoil ASA.
+
+  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 <http://www.gnu.org/licenses/>.
+*/
+
+#include "AutoDiff.hpp"
+
+#include <iostream>
+#include <cmath>
+
+struct Func
+{
+    template <typename T>
+    T operator()(T x) const
+    {
+#if 1
+        T r = std::sqrt(std::cos(x * x) + x) - 1.2;
+        return r;
+#else
+        return x;
+        // const int n = 6;
+        // double xv[6] = { 0.0, 0.2, 0.4, 0.6, 0.8, 1.0 };
+        // double yv[6] = { -0.5, -0.3, -0.1, 0.1, 0.3, 0.5 };
+        // int interv = -1;
+        // for (int i = 0; i < n; ++i) {
+        //     if (x < xv[i]) {
+        //         interv = i - 1;
+        //         break;
+        //     }
+        // }
+        // T t = (x - xv[interv])/(xv[interv+1] - xv[interv]);
+        // return (1.0 - t)*yv[interv] + t*yv[interv+1];
+#endif
+    }
+};
+
+// template <class ErrorPolicy = ThrowOnError>
+class Newton
+{
+public:
+    /// Implements a scalar Newton solve.
+    template <class Functor>
+    inline static double solve(const Functor& f,
+                               const double initial_guess,
+                               const int max_iter,
+                               const double tolerance,
+                               int& iterations_used)
+    {
+        double x = initial_guess;
+        iterations_used = 0;
+        typedef AutoDiff::Forward<double> AD;
+        while (std::abs(f(x)) > tolerance && ++iterations_used < max_iter) {
+            AD xfad = AD::variable(x);
+            AD rfad = f(xfad);
+            x = x - rfad.val()/rfad.der();
+        }
+        return x;
+    }
+};
+
+
+int main()
+{
+    int iter = 0;
+    const double atol = 1.0e-13;
+    const double soln = Newton::solve(Func(), 0.1, 30, atol, iter);
+
+    std::cout.precision(16);
+    std::cout << "Solution is: " << soln
+              << "   using " << iter << " iterations." << '\n';
+    std::cout << " f(x) = " << Func()(soln) << '\n';
+}

examples/sim_2p_incomp_adfi.cpp

@@ -0,0 +1,304 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
+
+
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif // HAVE_CONFIG_H
+
+#include <opm/core/pressure/FlowBCManager.hpp>
+
+#include <opm/core/grid.h>
+#include <opm/core/grid/GridManager.hpp>
+#include <opm/core/wells.h>
+#include <opm/core/wells/WellsManager.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/simulator/initState.hpp>
+#include <opm/core/simulator/SimulatorReport.hpp>
+#include <opm/core/simulator/SimulatorTimer.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+
+#include <opm/core/props/IncompPropertiesBasic.hpp>
+#include <opm/core/props/IncompPropertiesFromDeck.hpp>
+#include <opm/core/props/rock/RockCompressibility.hpp>
+
+#include <opm/core/linalg/LinearSolverFactory.hpp>
+
+#include <opm/core/simulator/TwophaseState.hpp>
+#include <opm/core/simulator/WellState.hpp>
+#include "SimulatorIncompTwophaseAdfi.hpp"
+
+#include <boost/scoped_ptr.hpp>
+#include <boost/filesystem.hpp>
+
+#include <algorithm>
+#include <iostream>
+#include <vector>
+#include <numeric>
+
+
+namespace
+{
+    void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
+    {
+        if (param.anyUnused()) {
+            std::cout << "--------------------   Unused parameters:   --------------------\n";
+            param.displayUsage();
+            std::cout << "----------------------------------------------------------------" << std::endl;
+        }
+    }
+} // anon namespace
+
+
+
+// ----------------- Main program -----------------
+int
+main(int argc, char** argv)
+{
+    using namespace Opm;
+
+    std::cout << "\n================    Test program for incompressible two-phase flow     ===============\n\n";
+    parameter::ParameterGroup param(argc, argv, false);
+    std::cout << "---------------    Reading parameters     ---------------" << std::endl;
+
+#if ! HAVE_SUITESPARSE_UMFPACK_H
+    // This is an extra check to intercept a potentially invalid request for the
+    // implicit transport solver as early as possible for the user.
+    {
+        const std::string transport_solver_type
+            = param.getDefault<std::string>("transport_solver_type", "ad");
+        if (transport_solver_type == "implicit") {
+            THROW("Cannot use implicit transport solver without UMFPACK. "
+                  "Either reconfigure opm-core with SuiteSparse/UMFPACK support and recompile, "
+                  "or use the reordering solver (transport_solver_type=reorder).");
+        }
+    }
+#endif
+
+    // If we have a "deck_filename", grid and props will be read from that.
+    bool use_deck = param.has("deck_filename");
+    boost::scoped_ptr<EclipseGridParser> deck;
+    boost::scoped_ptr<GridManager> grid;
+    boost::scoped_ptr<IncompPropertiesInterface> props;
+    boost::scoped_ptr<RockCompressibility> rock_comp;
+    TwophaseState state;
+    // bool check_well_controls = false;
+    // int max_well_control_iterations = 0;
+    double gravity[3] = { 0.0 };
+    if (use_deck) {
+        std::string deck_filename = param.get<std::string>("deck_filename");
+        deck.reset(new EclipseGridParser(deck_filename));
+        // Grid init
+        grid.reset(new GridManager(*deck));
+        // Rock and fluid init
+        props.reset(new IncompPropertiesFromDeck(*deck, *grid->c_grid()));
+        // check_well_controls = param.getDefault("check_well_controls", false);
+        // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
+        // Rock compressibility.
+        rock_comp.reset(new RockCompressibility(*deck));
+        // Gravity.
+        gravity[2] = deck->hasField("NOGRAV") ? 0.0 : unit::gravity;
+        // Init state variables (saturation and pressure).
+        if (param.has("init_saturation")) {
+            initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
+        } else {
+            initStateFromDeck(*grid->c_grid(), *props, *deck, gravity[2], state);
+        }
+    } else {
+        // Grid init.
+        const int nx = param.getDefault("nx", 100);
+        const int ny = param.getDefault("ny", 100);
+        const int nz = param.getDefault("nz", 1);
+        const double dx = param.getDefault("dx", 1.0);
+        const double dy = param.getDefault("dy", 1.0);
+        const double dz = param.getDefault("dz", 1.0);
+        grid.reset(new GridManager(nx, ny, nz, dx, dy, dz));
+        // Rock and fluid init.
+        props.reset(new IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
+        // Rock compressibility.
+        rock_comp.reset(new RockCompressibility(param));
+        // Gravity.
+        gravity[2] = param.getDefault("gravity", 0.0);
+        // Init state variables (saturation and pressure).
+        initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
+    }
+
+    // Warn if gravity but no density difference.
+    bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
+    if (use_gravity) {
+        if (props->density()[0] == props->density()[1]) {
+            std::cout << "**** Warning: nonzero gravity, but zero density difference." << std::endl;
+        }
+    }
+    const double *grav = use_gravity ? &gravity[0] : 0;
+
+    // Initialising src
+    int num_cells = grid->c_grid()->number_of_cells;
+    std::vector<double> src(num_cells, 0.0);
+    if (use_deck) {
+        // Do nothing, wells will be the driving force, not source terms.
+    } else {
+        // Compute pore volumes, in order to enable specifying injection rate
+        // terms of total pore volume.
+        std::vector<double> porevol;
+        if (rock_comp->isActive()) {
+            computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
+        } else {
+            computePorevolume(*grid->c_grid(), props->porosity(), porevol);
+        }
+        const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
+        const double default_injection = use_gravity ? 0.0 : 0.1;
+        const double flow_per_sec = param.getDefault<double>("injected_porevolumes_per_day", default_injection)
+            *tot_porevol_init/unit::day;
+        src[0] = flow_per_sec;
+        src[num_cells - 1] = -flow_per_sec;
+    }
+
+    // Boundary conditions.
+    FlowBCManager bcs;
+    if (param.getDefault("use_pside", false)) {
+        int pside = param.get<int>("pside");
+        double pside_pressure = param.get<double>("pside_pressure");
+        bcs.pressureSide(*grid->c_grid(), FlowBCManager::Side(pside), pside_pressure);
+    }
+
+    // Linear solver.
+    LinearSolverFactory linsolver(param);
+
+    // Write parameters used for later reference.
+    bool output = param.getDefault("output", true);
+    std::ofstream epoch_os;
+    std::string output_dir;
+    if (output) {
+        output_dir =
+            param.getDefault("output_dir", std::string("output"));
+        boost::filesystem::path fpath(output_dir);
+        try {
+            create_directories(fpath);
+        }
+        catch (...) {
+            THROW("Creating directories failed: " << fpath);
+        }
+        std::string filename = output_dir + "/epoch_timing.param";
+        epoch_os.open(filename.c_str(), std::fstream::trunc | std::fstream::out);
+        // open file to clean it. The file is appended to in SimulatorTwophase
+        filename = output_dir + "/step_timing.param";
+        std::fstream step_os(filename.c_str(), std::fstream::trunc | std::fstream::out);
+        step_os.close();
+        param.writeParam(output_dir + "/simulation.param");
+    }
+
+
+    std::cout << "\n\n================    Starting main simulation loop     ===============\n"
+              << "                        (number of epochs: "
+              << (use_deck ? deck->numberOfEpochs() : 1) << ")\n\n" << std::flush;
+
+    SimulatorReport rep;
+    if (!use_deck) {
+        // Simple simulation without a deck.
+        WellsManager wells; // no wells.
+        SimulatorIncompTwophaseAdfi simulator(param,
+                                          *grid->c_grid(),
+                                          *props,
+                                          rock_comp->isActive() ? rock_comp.get() : 0,
+                                          wells,
+                                          src,
+                                          bcs.c_bcs(),
+                                          linsolver,
+                                          grav);
+        SimulatorTimer simtimer;
+        simtimer.init(param);
+        warnIfUnusedParams(param);
+        WellState well_state;
+        well_state.init(0, state);
+        rep = simulator.run(simtimer, state, well_state);
+    } else {
+        // With a deck, we may have more epochs etc.
+        WellState well_state;
+        int step = 0;
+        SimulatorTimer simtimer;
+        // Use timer for last epoch to obtain total time.
+        deck->setCurrentEpoch(deck->numberOfEpochs() - 1);
+        simtimer.init(*deck);
+        const double total_time = simtimer.totalTime();
+        for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) {
+            // Set epoch index.
+            deck->setCurrentEpoch(epoch);
+
+            // Update the timer.
+            if (deck->hasField("TSTEP")) {
+                simtimer.init(*deck);
+            } else {
+                if (epoch != 0) {
+                    THROW("No TSTEP in deck for epoch " << epoch);
+                }
+                simtimer.init(param);
+            }
+            simtimer.setCurrentStepNum(step);
+            simtimer.setTotalTime(total_time);
+
+            // Report on start of epoch.
+            std::cout << "\n\n--------------    Starting epoch " << epoch << "    --------------"
+                      << "\n                  (number of steps: "
+                      << simtimer.numSteps() - step << ")\n\n" << std::flush;
+
+            // Create new wells, well_state
+            WellsManager wells(*deck, *grid->c_grid(), props->permeability());
+            // @@@ HACK: we should really make a new well state and
+            // properly transfer old well state to it every epoch,
+            // since number of wells may change etc.
+            if (epoch == 0) {
+                well_state.init(wells.c_wells(), state);
+            }
+
+            // Create and run simulator.
+            SimulatorIncompTwophaseAdfi simulator(param,
+                                              *grid->c_grid(),
+                                              *props,
+                                              rock_comp->isActive() ? rock_comp.get() : 0,
+                                              wells,
+                                              src,
+                                              bcs.c_bcs(),
+                                              linsolver,
+                                              grav);
+            if (epoch == 0) {
+                warnIfUnusedParams(param);
+            }
+            SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state);
+            if (output) {
+                epoch_rep.reportParam(epoch_os);
+            }
+            // Update total timing report and remember step number.
+            rep += epoch_rep;
+            step = simtimer.currentStepNum();
+        }
+    }
+
+    std::cout << "\n\n================    End of simulation     ===============\n\n";
+    rep.report(std::cout);
+
+    if (output) {
+      std::string filename = output_dir + "/walltime.param";
+      std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out);
+      rep.reportParam(tot_os);
+    }
+
+}

examples/sim_simple.cpp

@@ -0,0 +1,285 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
+
+#include "AutoDiffBlock.hpp"
+#include "AutoDiffHelpers.hpp"
+#include <opm/core/grid.h>
+#include <opm/core/grid/GridManager.hpp>
+#include <opm/core/props/IncompPropertiesBasic.hpp>
+#include <opm/core/utility/Units.hpp>
+#include <opm/core/utility/StopWatch.hpp>
+#include <opm/core/pressure/tpfa/trans_tpfa.h>
+#include <Eigen/UmfPackSupport>
+
+#include <iostream>
+#include <cstdlib>
+
+/*
+  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
+
+*/
+
+
+template <class ADB>
+std::vector<ADB>
+phaseMobility(const Opm::IncompPropertiesInterface& props,
+              const std::vector<int>& cells,
+              const typename ADB::V& sw)
+{
+    typedef Eigen::Array<double, Eigen::Dynamic, 2, Eigen::RowMajor> TwoCol;
+    typedef Eigen::Array<double, Eigen::Dynamic, 4, Eigen::RowMajor> 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();
+    std::vector<M> dmw = { krwjac/mu[0] };
+    std::vector<M> dmo = { krojac/mu[1] };
+
+    std::vector<ADB> pmobc = { ADB::function(krw / mu[0], dmw) ,
+                               ADB::function(kro / mu[1], dmo) };
+    return pmobc;
+}
+
+/// Returns fw(sw).
+template <class ADB>
+ADB
+fluxFunc(const std::vector<ADB>& m)
+{
+    assert (m.size() == 2);
+
+    ADB f = m[0] / (m[0] + m[1]);
+
+    return f;
+}
+
+
+int main()
+{
+    typedef AutoDiff::ForwardBlock<double> ADB;
+    typedef ADB::V V;
+    typedef ADB::M M;
+
+    Opm::time::StopWatch clock;
+    clock.start();
+    const Opm::GridManager gm(3,3);//(50, 50, 10);
+    const UnstructuredGrid& grid = *gm.c_grid();
+    using namespace Opm::unit;
+    using namespace Opm::prefix;
+    // 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);
+    V htrans(grid.cell_facepos[grid.number_of_cells]);
+    tpfa_htrans_compute(const_cast<UnstructuredGrid*>(&grid), props.permeability(), htrans.data());
+    V trans_all(grid.number_of_faces);
+    // tpfa_trans_compute(const_cast<UnstructuredGrid*>(&grid), htrans.data(), trans_all.data());
+    const int nc = grid.number_of_cells;
+    std::vector<int> allcells(nc);
+    for (int i = 0; i < nc; ++i) {
+        allcells[i] = i;
+    }
+    std::cerr << "Opm core " << clock.secsSinceLast() << std::endl;
+
+    // Define neighbourhood-derived operator matrices.
+    const HelperOps ops(grid);
+    const int num_internal = ops.internal_faces.size();
+    std::cerr << "Topology matrices " << clock.secsSinceLast() << std::endl;
+
+    typedef AutoDiff::ForwardBlock<double> ADB;
+    typedef ADB::V V;
+
+    // q
+    V q(nc);
+    q.setZero();
+    q[0] = 1.0;
+    q[nc-1] = -1.0;
+
+    // s0 - this is explicit now
+    typedef Eigen::Array<double, Eigen::Dynamic, 2, Eigen::RowMajor> TwoCol;
+    TwoCol s0(nc, 2);
+    s0.leftCols<1>().setZero();
+    s0.rightCols<1>().setOnes();
+
+    // totmob - explicit as well
+    TwoCol kr(nc, 2);
+    props.relperm(nc, s0.data(), allcells.data(), kr.data(), 0);
+    const V krw = kr.leftCols<1>();
+    const V kro = kr.rightCols<1>();
+    const double* mu = props.viscosity();
+    const V totmob = krw/mu[0] + kro/mu[1];
+
+    // Moved down here because we need total mobility.
+    tpfa_eff_trans_compute(const_cast<UnstructuredGrid*>(&grid), totmob.data(),
+                           htrans.data(), trans_all.data());
+    // Still explicit, and no upwinding!
+    V mobtransf(num_internal);
+    for (int fi = 0; fi < num_internal; ++fi) {
+        mobtransf[fi] = trans_all[ops.internal_faces[fi]];
+    }
+    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 variable.
+    const std::vector<int> bpat = { nc };
+    // Could actually write { nc } instead of bpat below,
+    // but we prefer a named variable since we will repeat it.
+    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.
+    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,
+    // 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<M> solver;
+    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;
+    }
+    // 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;
+    }
+    const V p1 = p0 - dp;
+    std::cerr << "Solve " << clock.secsSinceLast() << 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;
+    const V sw0 = s0.leftCols<1>();
+    // V sw1 = sw0;
+    V sw1 = 0.5*V::Ones(nc,1);
+    const V ndp = (ops.ngrad * p1.matrix()).array();
+    const V dflux = mobtransf * ndp;
+    const UpwindSelector<double> upwind(grid, ops, dflux);
+    const V pv = Eigen::Map<const V>(props.porosity(), nc, 1)
+        * Eigen::Map<const V>(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 ADB sw = ADB::variable(0, sw1, bpat);
+        const std::vector<ADB> pmobc = phaseMobility<ADB>(props, allcells, sw.value());
+        const std::vector<ADB> pmobf = upwind.select(pmobc);
+        const ADB fw_cell = fluxFunc(pmobc);
+        const ADB fw_face = fluxFunc(pmobf);
+        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;
+
+        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;
+        }
+        const V ds = solver.solve(transport_residual.value().matrix()).array();
+        if (solver.info() != Eigen::Success) {
+            std::cerr << "Transport solve failure\n";
+            return EXIT_FAILURE;
+        }
+        sw1 = sw.value() - ds;
+        std::cerr << "Solve for s[" << it << "]: "
+                  << clock.secsSinceLast() << '\n';
+        sw1 = sw1.min(V::Ones(nc,1)).max(V::Zero(nc,1));
+
+        it += 1;
+    } while (res_norm > 1e-7);
+
+    std::cout << "Saturation solution:\n"
+              << "function s1 = solution\n"
+              << "s1 = [\n" << sw1 << "\n];\n";
+}