C++ time::StopWatch类代码示例

// ----------------- Main program -----------------
int
main(int argc, char** argv)
try
{
    using namespace Opm;
    if (argc <= 1) {
        usage();
        exit(1);
    } 
    const char* eclipseFilename = argv[1];
    EclipseStateConstPtr eclState; 
    ParserPtr parser(new Opm::Parser);
    Opm::ParseContext parseContext({{ ParseContext::PARSE_RANDOM_SLASH , InputError::IGNORE }, 
                              { ParseContext::PARSE_UNKNOWN_KEYWORD, InputError::IGNORE},
                              { ParseContext::PARSE_RANDOM_TEXT, InputError::IGNORE},
                              { ParseContext::UNSUPPORTED_SCHEDULE_GEO_MODIFIER, InputError::IGNORE},
                              { ParseContext::UNSUPPORTED_COMPORD_TYPE, InputError::IGNORE},
                              { ParseContext::UNSUPPORTED_INITIAL_THPRES, InputError::IGNORE},
                              { ParseContext::INTERNAL_ERROR_UNINITIALIZED_THPRES, InputError::IGNORE}
                             });
    Opm::DeckConstPtr deck(parser->parseFile(eclipseFilename, parseContext));
    eclState.reset(new EclipseState(deck, parseContext));

    GridManager gm(deck);
    const UnstructuredGrid& grid = *gm.c_grid();
    using boost::filesystem::path; 
    path fpath(eclipseFilename);
    std::string baseName;
    if (boost::to_upper_copy(path(fpath.extension()).string())== ".DATA") {
        baseName = path(fpath.stem()).string();
    } else {
        baseName = path(fpath.filename()).string();
    }

    std::string logFile = baseName + ".SATFUNCLOG";
    Opm::time::StopWatch timer;
    timer.start();
    RelpermDiagnostics diagnostic(logFile);
    diagnostic.diagnosis(eclState, deck, grid);
    timer.stop();
    double tt = timer.secsSinceStart();
    std::cout << "relperm diagnostics: " << tt << " seconds." << std::endl;
}
catch (const std::exception &e) {
    std::cerr << "Program threw an exception: " << e.what() << "\n";
    throw;
}

void test_flowsolver(const GI& g, const RI& r, double tol, int kind)
{
    typedef typename GI::CellIterator                   CI;
    typedef typename CI::FaceIterator                   FI;
    typedef double (*SolutionFuncPtr)(const Vec&);

    //typedef Opm::BasicBoundaryConditions<true, false>  FBC;
    typedef Opm::FunctionBoundaryConditions<SolutionFuncPtr> FBC;
    typedef Opm::IncompFlowSolverHybrid<GI, RI, FBC,
        Opm::MimeticIPEvaluator> FlowSolver;

    FlowSolver solver;

    // FBC flow_bc;
    // assign_bc(g, flow_bc);
    FBC flow_bc(&u);

    typename CI::Vector gravity(0.0);

    std::cout << "========== Init pressure solver =============" << std::endl;
    Opm::time::StopWatch rolex;
    rolex.start();
    solver.init(g, r, gravity, flow_bc);
    rolex.stop();
    std::cout << "========== Time in seconds: " << rolex.secsSinceStart() << " =============" << std::endl;

    std::vector<double> src(g.numberOfCells(), 0.0);
    assign_src(g, src);
    std::vector<double> sat(g.numberOfCells(), 0.0);


    std::cout << "========== Starting pressure solve =============" << std::endl;
    rolex.start();
    solver.solve(r, sat, flow_bc, src, tol, 3, kind);
    rolex.stop();
    std::cout << "========== Time in seconds: " << rolex.secsSinceStart() << " =============" << std::endl;

    typedef typename FlowSolver::SolutionType FlowSolution;
    FlowSolution soln = solver.getSolution();

    std::vector<typename GI::Vector> cell_velocity;
    estimateCellVelocity(cell_velocity, g, solver.getSolution());
    // Dune's vtk writer wants multi-component data to be flattened.
    std::vector<double> cell_velocity_flat(&*cell_velocity.front().begin(),
                                           &*cell_velocity.back().end());
    std::vector<double> cell_pressure;
    getCellPressure(cell_pressure, g, soln);

    compare_pressure(g, cell_pressure);

    Dune::VTKWriter<typename GI::GridType::LeafGridView> vtkwriter(g.grid().leafView());
    vtkwriter.addCellData(cell_velocity_flat, "velocity", GI::GridType::dimension);
    vtkwriter.addCellData(cell_pressure, "pressure");
    vtkwriter.write("testsolution-" + boost::lexical_cast<std::string>(0),
                    Dune::VTKOptions::ascii);
}

//! \brief Main driver
int main(int argc, char** argv)
try
{
    try {
        static const int dim = 3;

        typedef Dune::CpGrid GridType;

        if (argc < 2 || strcmp(argv[1],"-h") == 0
                || strcmp(argv[1],"--help") == 0
                || strcmp(argv[1],"-?") == 0) {
            syntax(argv);
            exit(1);
        }
        Params p;
        parseCommandLine(argc,argv,p);

        Opm::time::StopWatch watch;
        watch.start();

        GridType grid;
        if (p.file == "uniform") {
            std::array<int,3> cells;
            cells[0] = p.cellsx;
            cells[1] = p.cellsy;
            cells[2] = p.cellsz;
            std::array<double,3> cellsize;
            cellsize[0] = cellsize[1] = cellsize[2] = 1.f;
            grid.createCartesian(cells,cellsize);
        } else
            grid.readEclipseFormat(p.file,p.ctol,false);

        typedef GridType::ctype ctype;
        Opm::Elasticity::ElasticityUpscale<GridType> upscale(grid, p.ctol,
                p.Emin, p.file,
                p.rocklist,
                p.verbose);
        if (p.max[0] < 0 || p.min[0] < 0) {
            std::cout << "determine side coordinates..." << std::endl;
            upscale.findBoundaries(p.min,p.max);
            std::cout << "  min " << p.min[0] << " " << p.min[1] << " " << p.min[2] << std::endl;
            std::cout << "  max " << p.max[0] << " " << p.max[1] << " " << p.max[2] << std::endl;
        }
        if (p.n1 == -1 || p.n2 == -1) {
            p.n1 = grid.logicalCartesianSize()[0];
            p.n2 = grid.logicalCartesianSize()[1];
        }
        if (p.linsolver.zcells == -1) {
            double lz = p.max[2]-p.min[2];
            int nz = grid.logicalCartesianSize()[2];
            double hz = lz/nz;
            double lp = sqrt((double)(p.max[0]-p.min[0])*(p.max[1]-p.min[1]));
            int np = std::max(grid.logicalCartesianSize()[0],
                              grid.logicalCartesianSize()[1]);
            double hp = lp/np;
            p.linsolver.zcells = (int)(2*hp/hz+0.5);
        }
        std::cout << "logical dimension: " << grid.logicalCartesianSize()[0]
                  << "x"                   << grid.logicalCartesianSize()[1]
                  << "x"                   << grid.logicalCartesianSize()[2]
                  << std::endl;

        if (p.method == UPSCALE_MPC) {
            std::cout << "using MPC couplings in all directions..." << std::endl;
            upscale.periodicBCs(p.min, p.max);
            std::cout << "preprocessing grid..." << std::endl;
            upscale.A.initForAssembly();
        } else if (p.method == UPSCALE_MORTAR) {
            std::cout << "using Mortar couplings.." << std::endl;
            upscale.periodicBCsMortar(p.min, p.max, p.n1, p.n2,
                                      p.lambda[0], p.lambda[1]);
        } else if (p.method == UPSCALE_NONE) {
            std::cout << "no periodicity approach applied.." << std::endl;
            upscale.fixCorners(p.min, p.max);
            upscale.A.initForAssembly();
        }
        Dune::FieldMatrix<double,6,6> C;
        Dune::VTKWriter<GridType::LeafGridView>* vtkwriter=0;
        if (!p.vtufile.empty())
            vtkwriter = new Dune::VTKWriter<GridType::LeafGridView>(grid.leafView());
        Opm::Elasticity::Vector field[6];
        std::cout << "assembling elasticity operator..." << "\n";
        upscale.assemble(-1,true);
        std::cout << "setting up linear solver..." << std::endl;
        upscale.setupSolvers(p.linsolver);

//#pragma omp parallel for schedule(static)
        for (int i=0; i<6; ++i) {
            std::cout << "processing case " << i+1 << "..." << std::endl;
            std::cout << "\tassembling load vector..." << std::endl;
            upscale.assemble(i,false);
            std::cout << "\tsolving..." << std::endl;
            upscale.solve(i);
            upscale.A.expandSolution(field[i],upscale.u[i]);
#define CLAMP(x) (fabs(x)<1.e-4?0.0:x)
            for (size_t j=0; j<field[i].size(); ++j) {
                double val = field[i][j];
                field[i][j] = CLAMP(val);
            }
//.........这里部分代码省略.........

//! \brief Write a log of the simulation to a text file
void writeOutput(const Params& p, Opm::time::StopWatch& watch, int cells,
                 const std::vector<double>& volume,
                 const Dune::FieldMatrix<double,6,6>& C)
{
    // get current time
    time_t rawtime;
    struct tm* timeinfo;
    time(&rawtime);
    timeinfo = localtime(&rawtime);

    // get hostname
    char hostname[1024];
    gethostname(hostname,1024);

    std::string method = "mortar";
    if (p.method == UPSCALE_MPC)
        method = "mpc";
    if (p.method == UPSCALE_NONE)
        method = "none";

    // write log
    std::ofstream f;
    f.open(p.output.c_str());
    f << "######################################################################" << std::endl
      << "# Results from upscaling elastic moduli." << std::endl
      << "#" << std::endl
      << "# Finished: " << asctime(timeinfo)
      << "# Hostname: " << hostname << std::endl
      << "#" << std::endl
      << "# Upscaling time: " << watch.secsSinceStart() << " secs" << std::endl
      << "#" << std::endl;
    if (p.file == "uniform") {
        f << "# Uniform grid used" << std::endl
          << "#\t cells: " << p.cellsx*p.cellsy*p.cellsz << std::endl;
    }
    else {
        f  << "# Eclipse file: " << p.file << std::endl
           << "#\t cells: " << cells << std::endl;
    }
    f << "#" << std::endl;
    if (!p.rocklist.empty()) {
        f << "# Rock list: " << p.rocklist << std::endl
          << "#" << std::endl;
    }
    f << "# Options used:" << std::endl
      << "#\t         method: " << method << std::endl
      << "#\t linsolver_type: " << (p.linsolver.type==Opm::Elasticity::DIRECT?"direct":"iterative")
      << std::endl;
    if (p.linsolver.type == Opm::Elasticity::ITERATIVE)
        f << "#\t           ltol: " << p.linsolver.tol << std::endl;
    if (p.file == "uniform") {
        f << "#\t          cellsx: " << p.cellsx << std::endl
          << "#\t          cellsy: " << p.cellsy << std::endl
          << "#\t          cellsz: " << p.cellsz << std::endl;
    }
    f << "#" << std::endl
      <<"# Materials: " << volume.size() << std::endl;
    for (size_t i=0; i<volume.size(); ++i)
        f << "#\t Material" << i+1 << ": " << volume[i]*100 << "%" << std::endl;
    f << "#" << std::endl
      << "######################################################################" << std::endl
      << C << std::endl;
}

    bool EulerUpstreamImplicit<GI, RP, BC>::transportSolve(std::vector<double>& saturation,
                                                           const double time,
                                                           const typename GI::Vector& /* gravity */,
                                                           const PressureSolution& pressure_sol,
                                                           const Opm::SparseVector<double>& /* injection_rates */) const
    {

        Opm::ReservoirState<2> state(mygrid_.c_grid());
        {
            std::vector<double>& sat = state.saturation();
            for (int i=0; i < mygrid_.numCells(); ++i){
                sat[2*i] = saturation[i];
                sat[2*i+1] = 1-saturation[i];
            }
        }

        //int count=0;
        const UnstructuredGrid* cgrid = mygrid_.c_grid();
        int numhf = cgrid->cell_facepos[cgrid->number_of_cells];

        std::vector<double>     faceflux(numhf);

        for (int c = 0, i = 0; c < cgrid->number_of_cells; ++c){
            for (; i < cgrid->cell_facepos[c + 1]; ++i) {
                int f= cgrid->cell_faces[i];
                double outflux = pressure_sol.outflux(i);
                double sgn = 2.0*(cgrid->face_cells[2*f + 0] == c) - 1;
                faceflux[f] = sgn * outflux;
            }
        }
        int num_db=direclet_hfaces_.size();
        std::vector<double> sflux(num_db);
        for (int i=0; i < num_db;++i){
            sflux[i]=-pressure_sol.outflux(direclet_hfaces_[i]);
        }
        state.faceflux()=faceflux;

        double dt_transport = time;
        int nr_transport_steps = 1;
	Opm::time::StopWatch clock;
        int repeats = 0;
        bool finished = false;
        clock.start();

        TwophaseFluid myfluid(myrp_);
        double* tmp_grav=0;
        const UnstructuredGrid& c_grid=*mygrid_.c_grid();
        TransportModel model(myfluid,c_grid,porevol_,tmp_grav);
        model.makefhfQPeriodic(periodic_faces_,periodic_hfaces_, periodic_nbfaces_);
        model.initGravityTrans(*mygrid_.c_grid(),htrans_);
        TransportSolver tsolver(model);
        LinearSolver linsolve_;
        Opm::ImplicitTransportDetails::NRReport  rpt_;

        Opm::TransportSource tsrc;//create_transport_source(0, 2);
        // the input flux is assumed to be the saturation times the flux in the transport solver

        tsrc.nsrc =direclet_cells_.size();
        tsrc.saturation = direclet_sat_;
        tsrc.cell = direclet_cells_;
        tsrc.flux = sflux;

        while (!finished) {
            for (int q = 0; q < nr_transport_steps; ++q) {
                tsolver.solve(*mygrid_.c_grid(), &tsrc, dt_transport, ctrl_, state, linsolve_, rpt_);
                if(rpt_.flag<0){
                    break;
                }
            }
            if(!(rpt_.flag<0) ){
                finished =true;
            }else{
                if(repeats >max_repeats_){
                    finished=true;
                }else{
                    OPM_MESSAGE("Warning: Transport failed, retrying with more steps.");
                    nr_transport_steps *= 2;
                    dt_transport = time/nr_transport_steps;
                    if (ctrl_.verbosity){
                        std::cout << "Warning: Transport failed, retrying with more steps. dt = "
                                  << dt_transport/Opm::unit::year << " year.\n";
                    }

                    std::vector<double>& sat = state.saturation();
                    for (int i=0; i < mygrid_.numCells(); ++i){
                        sat[2*i] = saturation[i];
                        sat[2*i+1] = 1-saturation[i];
                    }
                }
            }
            repeats +=1;
        }
        clock.stop();
        std::cout << "EulerUpstreamImplicite used  " << repeats
                  << " repeats and " << nr_transport_steps <<" steps"<< std::endl;
#ifdef VERBOSE
        std::cout << "Seconds taken by transport solver: " << clock.secsSinceStart() << std::endl;
#endif // VERBOSE
        {
            std::vector<double>& sat = state.saturation();
//.........这里部分代码省略.........

CollOfScalar EquelleRuntimeCUDA::newtonSolve(const ResidualFunctor& rescomp,
                                            const CollOfScalar& u_initialguess)
{
    Opm::time::StopWatch clock;
    clock.start();

    // Set up Newton loop.
 
    // Define the primary variable
    CollOfScalar u = CollOfScalar(u_initialguess, true);
 
    if (verbose_ > 2) {
        output("Initial u", u);
        output("    newtonSolve: norm (initial u)", twoNorm(u));
    }
    CollOfScalar residual = rescomp(u);   
    if (verbose_ > 2) {
        output("Initial residual", residual);
        output("    newtonSolve: norm (initial residual)", twoNorm(residual));
    }

    int iter = 0;

    // Debugging output not specified in Equelle.
    if (verbose_ > 1) {
        std::cout << "    newtonSolve: iter = " << iter << " (max = " << max_iter_
		  << "), norm(residual) = " << twoNorm(residual)
                  << " (tol = " << abs_res_tol_ << ")" << std::endl;
    }

    CollOfScalar du;

    // Execute newton loop until residual is small or we have used too many iterations.
    while ( (twoNorm(residual) > abs_res_tol_) && (iter < max_iter_) ) {
	
	if ( solver_.getSolver() == CPU ) {
	    du = serialSolveForUpdate(residual);
	}
	else {
	    // Solve linear equations for du, apply update.
	    du = solver_.solve(residual.derivative(),
			       residual.value(),
			       verbose_);
	}

	// du is a constant, hence, u is still a primary variable with an identity
	// matrix as its derivative.
	u = u - du;

        // Recompute residual.
        residual = rescomp(u);

        if (verbose_ > 2) {
            // Debugging output not specified in Equelle.
            output("u", u);
            output("    newtonSolve: norm(u)", twoNorm(u));
            output("residual", residual);
            output("    newtonSolve: norm(residual)", twoNorm(residual));
        }

        ++iter;

        // Debugging output not specified in Equelle.
        if (verbose_ > 1) {
            std::cout << "    newtonSolve: iter = " << iter << " (max = " << max_iter_
		      << "), norm(residual) = " << twoNorm(residual)
                      << " (tol = " << abs_res_tol_ << ")" << std::endl;
        }

    }
    if (verbose_ > 0) {
        if (twoNorm(residual) > abs_res_tol_) {
            std::cout << "Newton solver failed to converge in " << max_iter_ << " iterations" << std::endl;
        } else {
            std::cout << "Newton solver converged in " << iter << " iterations" << std::endl;
        }
    }

    if (verbose_ > 1) {
        std::cout << "Newton solver took: " << clock.secsSinceLast() << " seconds." << std::endl;
    }

    return CollOfScalar(u.value());
}

        void buildFaceIndices()
        {
#ifdef VERBOSE
            std::cout << "Building unique face indices... " << std::flush;
	    Opm::time::StopWatch clock;
            clock.start();
#endif
            typedef CellIterator CI;
            typedef typename CI::FaceIterator FI;

            // We build the actual cell to face mapping in two passes.
            // [code mostly lifted from IncompFlowSolverHybrid::enumerateGridDof(),
            //  but with a twist: This code builds a mapping from cells in index
            //  order to unique face numbers, while the mapping built in the
            //  enumerateGridDof() method was ordered by cell iterator order]

            // Allocate and reserve structures.
            const int nc = numberOfCells();
            std::vector<int> cell(nc, -1);
            std::vector<int> num_faces(nc); // In index order.
            std::vector<int> fpos;     fpos  .reserve(nc + 1);
            std::vector<int> num_cf;   num_cf.reserve(nc); // In iterator order.
            std::vector<int> faces ;

            // First pass: enumerate internal faces.
            int cellno = 0; fpos.push_back(0);
            int tot_ncf = 0, tot_ncf2 = 0, max_ncf = 0;
            for (CI c = cellbegin(); c != cellend(); ++c, ++cellno) {
                const int c0 = c->index();
                ASSERT((0 <= c0) && (c0 < nc) && (cell[c0] == -1));
                cell[c0] = cellno;
                num_cf.push_back(0);
                int& ncf = num_cf.back();
                for (FI f = c->facebegin(); f != c-> faceend(); ++f) {
                    if (!f->boundary()) {
                        const int c1 = f->neighbourCellIndex();
                        ASSERT((0 <= c1) && (c1 < nc) && (c1 != c0));

                        if (cell[c1] == -1) {
                            // Previously undiscovered internal face.
                            faces.push_back(c1);
                        }
                    }
                    ++ncf;
                }
                num_faces[c0] = ncf;
                fpos.push_back(int(faces.size()));
                max_ncf  = std::max(max_ncf, ncf);
                tot_ncf  += ncf;
                tot_ncf2 += ncf * ncf;
            }
            ASSERT(cellno == nc);

            // Build cumulative face sizes enabling direct insertion of
            // face indices into cfdata later.
            std::vector<int> cumul_num_faces(numberOfCells() + 1);
            cumul_num_faces[0] = 0;
            std::partial_sum(num_faces.begin(), num_faces.end(), cumul_num_faces.begin() + 1);

            // Avoid (most) allocation(s) inside 'c' loop.
            std::vector<int>    l2g;
            l2g.reserve(max_ncf);
            std::vector<double> cfdata(tot_ncf);
            int total_num_faces = int(faces.size());

            // Second pass: build cell-to-face mapping, including boundary.
            typedef std::vector<int>::iterator VII;
            for (CI c = cellbegin(); c != cellend(); ++c) {
                const int c0 = c->index();
                ASSERT ((0 <=      c0 ) && (     c0  < nc) &&
                        (0 <= cell[c0]) && (cell[c0] < nc));
                const int ncf = num_cf[cell[c0]];
                l2g.resize(ncf, 0);
                for (FI f = c->facebegin(); f != c->faceend(); ++f) {
                    if (f->boundary()) {
                        // External, not counted before.  Add new face...
                        l2g[f->localIndex()] = total_num_faces++;
                    } else {
                        // Internal face.  Need to determine during
                        // traversal of which cell we discovered this
                        // face first, and extract the face number
                        // from the 'faces' table range of that cell.

                        // Note: std::find() below is potentially
                        // *VERY* expensive (e.g., large number of
                        // seeks in moderately sized data in case of
                        // faulted cells).
                        const int c1 = f->neighbourCellIndex();
                        ASSERT ((0 <=      c1 ) && (     c1  < nc) &&
                                (0 <= cell[c1]) && (cell[c1] < nc));

                        int t = c0, seek = c1;
                        if (cell[seek] < cell[t])
                            std::swap(t, seek);
                        int s = fpos[cell[t]], e = fpos[cell[t] + 1];
                        VII p = std::find(faces.begin() + s, faces.begin() + e, seek);
                        ASSERT(p != faces.begin() + e);
                        l2g[f->localIndex()] = p - faces.begin();
                    }
                }
//.........这里部分代码省略.........

    void EulerUpstream<GI, RP, BC>::transportSolve(std::vector<double>& saturation,
						   const double time,
						   const typename GI::Vector& gravity,
						   const PressureSolution& pressure_sol,
						   const Opm::SparseVector<double>& injection_rates) const
    {
	// Compute the cfl time-step.
	double cfl_dt = computeCflTime(saturation, time, gravity, pressure_sol);

	// Compute the number of small steps to take, and the actual small timestep.
	int nr_transport_steps;
	if (cfl_dt > time){
	    nr_transport_steps = minimum_small_steps_;
	} else {
            double steps = std::min<double>(std::ceil(time/cfl_dt), std::numeric_limits<int>::max());
            nr_transport_steps = int(steps);
	    nr_transport_steps = std::max(nr_transport_steps, minimum_small_steps_);
            nr_transport_steps = std::min(nr_transport_steps, maximum_small_steps_);
	}
	double dt_transport = time/nr_transport_steps;

	// Do the timestepping. The try-catch blocks are there to handle
	// the situation that smallTimeStep throws, which may happen due
	// to saturation out of bounds (if check_sat_ is true).
	// We cannot guarantee that this does not happen, since we do not
	// (yet) compute a capillary cfl condition.
	// Using exception for "alternate control flow" like this is bad
	// design, should rather use error return values for this.
	std::vector<double> saturation_initial(saturation);
	bool finished = false;
	int repeats = 0;
	const int max_repeats = 10;
	Opm::time::StopWatch clock;
        clock.start();
	while (!finished) {
	    try {
#ifdef VERBOSE
		std::cout << "Doing " << nr_transport_steps
			  << " steps for saturation equation with stepsize "
			  << dt_transport << " in seconds." << std::endl;
#endif // VERBOSE
		for (int q = 0; q < nr_transport_steps; ++q) {
		    smallTimeStep(saturation,
				  dt_transport,
				  gravity,
				  pressure_sol,
                                  injection_rates);
		}
		finished = true;
	    }
	    catch (...) {
		++repeats;
		if (repeats > max_repeats) {
		    throw;
		}
		OPM_MESSAGE("Warning: Transport failed, retrying with more steps.");
		nr_transport_steps *= 2;
		dt_transport = time/nr_transport_steps;
		saturation = saturation_initial;
	    }
	}
        clock.stop();
#ifdef VERBOSE
        std::cout << "Seconds taken by transport solver: " << clock.secsSinceStart() << std::endl;
#endif // VERBOSE
    }

    void ImplicitCapillarity<GI, RP, BC, IP>::transportSolve(std::vector<double>& saturation,
                                                             const double /*time*/,
                                                             const typename GI::Vector& gravity,
                                                             const PressureSolution& pressure_sol,
                                                             const Opm::SparseVector<double>& injection_rates) const
    {
        // Start a timer.
	Opm::time::StopWatch clock;
        clock.start();

        // Compute capillary mobilities.
        typedef typename RP::Mobility Mob;
        int num_cells = saturation.size();
        std::vector<Mob> cap_mob(num_cells);
        for (int c = 0; c < num_cells; ++c) {
            Mob& m = cap_mob[c];
            residual_.reservoirProperties().phaseMobility(0, c, saturation[c], m.mob);
            Mob mob2;
            residual_.reservoirProperties().phaseMobility(1, c, saturation[c], mob2.mob);
            Mob mob_tot;
            mob_tot.setToSum(m, mob2);
            Mob mob_totinv;
            mob_totinv.setToInverse(mob_tot);
            m *= mob_totinv;
            m *= mob2;
            ImplicitCapillarityDetails::thresholdMobility(m.mob, 1e-10); // @@TODO: User-set limit.
            // std::cout << m.mob(0,0) << '\n';
        }
        ReservoirPropertyFixedMobility<Mob> capillary_mobilities(cap_mob);

        // Set up boundary conditions.
        BC cap_press_bcs(residual_.boundaryConditions());
        for (int i = 0; i < cap_press_bcs.size(); ++i) {
            if (cap_press_bcs.flowCond(i).isPeriodic()) {
                cap_press_bcs.flowCond(i) = FlowBC(FlowBC::Periodic, 0.0);
            }
        }

        // Compute injection rates from residual.
        std::vector<double> injection_rates_residual(num_cells);
	residual_.computeResidual(saturation, gravity, pressure_sol, injection_rates,
                                  method_viscous_, method_gravity_, false,
                                  injection_rates_residual);
        for (int i = 0; i < num_cells; ++i) {
            injection_rates_residual[i] = -injection_rates_residual[i];
        }

        // Compute capillary pressure.
        // Note that the saturation is just a dummy for this call, since the mobilities are fixed.
        psolver_.solve(capillary_mobilities, saturation, cap_press_bcs, injection_rates_residual,
                       residual_tolerance_, linsolver_verbosity_, linsolver_type_);

        // Solve for constant to change capillary pressure solution by.
        std::vector<double> cap_press(num_cells);
        const PressureSolution& pcapsol = psolver_.getSolution();
        for (CIt c = residual_.grid().cellbegin(); c != residual_.grid().cellend(); ++c) {
            cap_press[c->index()] = pcapsol.pressure(c);
        }
        MatchSaturatedVolumeFunctor<GI, RP> functor(residual_.grid(),
                                                    residual_.reservoirProperties(),
                                                    saturation,
                                                    cap_press);
        double min_cap_press = *std::min_element(cap_press.begin(), cap_press.end());
        double max_cap_press = *std::max_element(cap_press.begin(), cap_press.end());
        double cap_press_range = max_cap_press - min_cap_press;
        double mod_low = 1e100;
        double mod_high = -1e100;
	Opm::bracketZero(functor, 0.0, cap_press_range, mod_low, mod_high);
        const int max_iter = 40;
        const double nonlinear_tolerance = 1e-12;
        int iterations_used = -1;
        typedef Opm::RegulaFalsi<Opm::ThrowOnError> RootFinder;
        double mod_correct = RootFinder::solve(functor, mod_low, mod_high, max_iter, nonlinear_tolerance, iterations_used);
        std::cout << "Moved capillary pressure solution by " << mod_correct << " after "
                  << iterations_used << " iterations." << std::endl;
        // saturation = functor.lastSaturations();
        const std::vector<double>& sat_new = functor.lastSaturations();
        for (int i = 0; i < num_cells; ++i) {
            saturation[i] = (1.0 - update_relaxation_)*saturation[i] + update_relaxation_*sat_new[i];
        }

        // Optionally check and/or clamp results.
	if (check_sat_ || clamp_sat_) {
	    checkAndPossiblyClampSat(saturation);
	}

        // Stop timer and optionally print seconds taken.
        clock.stop();
#ifdef VERBOSE
        std::cout << "Seconds taken by transport solver: " << clock.secsSinceStart() << std::endl;
#endif // VERBOSE
    }

int main()
try
{
    typedef Opm::AutoDiffBlock<double> ADB;
    typedef ADB::V V;
    typedef Eigen::SparseMatrix<double> S;

    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 Opm::HelperOps ops(grid);
    const int num_internal = ops.internal_faces.size();
    std::cerr << "Topology matrices " << clock.secsSinceLast() << std::endl;

    typedef Opm::AutoDiffBlock<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].

//.........这里部分代码省略.........

int run(Params& p)
{
  try {
    static const int dim = 3;

    Opm::time::StopWatch watch;
    watch.start();

    GridType grid;
    if (p.file == "uniform") {
      std::array<int,3> cells;
      cells[0] = p.cellsx;
      cells[1] = p.cellsy;
      cells[2] = p.cellsz;
      std::array<double,3> cellsize;
      cellsize[0] = p.max[0] > -1?p.max[0]/cells[0]:1.0;
      cellsize[1] = p.max[1] > -1?p.max[1]/cells[1]:1.0;
      cellsize[2] = p.max[2] > -1?p.max[2]/cells[2]:1.0;
      grid.createCartesian(cells,cellsize);
    } else {
        Opm::ParseContext parseContext;
        Opm::ParserPtr parser(new Opm::Parser());
        Opm::DeckConstPtr deck(parser->parseFile(p.file , parseContext));
        Opm::EclipseGrid inputGrid(deck);
        grid.processEclipseFormat(inputGrid, false);
    }
    ElasticityUpscale<GridType, AMG> upscale(grid, p.ctol, p.Emin, p.file,
                                             p.rocklist, p.verbose);

    if (p.max[0] < 0 || p.min[0] < 0) {
      std::cout << "determine side coordinates..." << std::endl;
      upscale.findBoundaries(p.min,p.max);
      std::cout << "  min " << p.min[0] << " " << p.min[1] << " " << p.min[2] << std::endl;
      std::cout << "  max " << p.max[0] << " " << p.max[1] << " " << p.max[2] << std::endl;
    }
    if (p.n1 == -1 || p.n2 == -1) {
      p.n1 = grid.logicalCartesianSize()[0];
      p.n2 = grid.logicalCartesianSize()[1];
    }
    if (p.linsolver.zcells == -1) {
      double lz = p.max[2]-p.min[2];
      int nz = grid.logicalCartesianSize()[2];
      double hz = lz/nz;
      double lp = sqrt((double)(p.max[0]-p.min[0])*(p.max[1]-p.min[1]));
      int np = std::max(grid.logicalCartesianSize()[0],
                        grid.logicalCartesianSize()[1]);
      double hp = lp/np;
      p.linsolver.zcells = (int)(2*hp/hz+0.5);
    }
    std::cout << "logical dimension: " << grid.logicalCartesianSize()[0]
              << "x"                   << grid.logicalCartesianSize()[1]
              << "x"                   << grid.logicalCartesianSize()[2]
              << std::endl;

    if (p.inspect == "mesh")
      return 0;

    if (p.linsolver.pre == UNDETERMINED) {
      double hx = (p.max[0]-p.min[0])/grid.logicalCartesianSize()[0];
      double hy = (p.max[1]-p.min[1])/grid.logicalCartesianSize()[1];
      double hz = (p.max[2]-p.min[2])/grid.logicalCartesianSize()[2];
      double aspect = sqrt(hx*hy)/hz;
      std::cout << "Estimated cell aspect ratio: " << aspect;
      if (aspect > 80) {
        p.linsolver.pre = TWOLEVEL;
        std::cout << " => Using two level preconditioner" << std::endl;
      } else {
        p.linsolver.pre = Opm::Elasticity::AMG;
        std::cout << " => Using AMG preconditioner" << std::endl;
      }
    }

    if (p.method == UPSCALE_MPC) {
      std::cout << "using MPC couplings in all directions..." << std::endl;
      upscale.periodicBCs(p.min, p.max);
      std::cout << "preprocessing grid..." << std::endl;
      upscale.A.initForAssembly();
    } else if (p.method == UPSCALE_MORTAR) {
      std::cout << "using Mortar couplings.." << std::endl;
      upscale.periodicBCsMortar(p.min, p.max, p.n1, p.n2,
                                p.lambda[0], p.lambda[1]);
    } else if (p.method == UPSCALE_NONE) {
      std::cout << "no periodicity approach applied.." << std::endl;
      upscale.fixCorners(p.min, p.max);
      upscale.A.initForAssembly();
    }

    Dune::FieldMatrix<double,6,6> C;
    Opm::Elasticity::Vector field[6];
    std::cout << "assembling elasticity operator..." << "\n";
    upscale.assemble(-1,true);
    std::cout << "setting up linear solver..." << std::endl;
    upscale.setupSolvers(p.linsolver);

    if (p.inspect == "load")
      Dune::storeMatrixMarket(upscale.A.getOperator(), "A.mtx");

    // the uzawa solver cannot run multithreaded
#ifdef HAVE_OPENMP
    if (p.linsolver.uzawa) {
//.........这里部分代码省略.........