C++ StopWatch::secsSinceStart方法代码示例

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);
}

// ----------------- 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;
}

//! \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();
//.........这里部分代码省略.........

//.........这里部分代码省略.........
            //  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();
                    }
                }
                ASSERT(int(l2g.size()) == num_faces[c0]);
                std::copy(l2g.begin(), l2g.end(), cfdata.begin() + cumul_num_faces[c0]);
            }
            num_faces_ = total_num_faces;
            max_faces_per_cell_ = max_ncf;
            face_indices_.assign(cfdata.begin(), cfdata.end(),
                                 num_faces.begin(), num_faces.end());

#ifdef VERBOSE
            clock.stop();
            double elapsed = clock.secsSinceStart();
            std::cout << "done.     Time elapsed: " << elapsed << std::endl;
#endif
        }

    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
    }