C++ SimulatorTimer::currentStepNum方法代码示例

    SimulatorReport SimulatorCompressiblePolymer::Impl::run(SimulatorTimer& timer,
                                                            PolymerBlackoilState& state,
                                                            WellState& well_state)
    {
        std::vector<double> transport_src(grid_.number_of_cells);
        std::vector<double> polymer_inflow_c(grid_.number_of_cells);

        // Initialisation.
        std::vector<double> initial_pressure;
        std::vector<double> porevol;
        if (rock_comp_props_ && rock_comp_props_->isActive()) {
            computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
        } else {
            computePorevolume(grid_, props_.porosity(), porevol);
        }
        const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
        std::vector<double> initial_porevol = porevol;

        // Main simulation loop.
        Opm::time::StopWatch pressure_timer;
        double ptime = 0.0;
        Opm::time::StopWatch transport_timer;
        double ttime = 0.0;
        Opm::time::StopWatch total_timer;
        total_timer.start();
        double init_surfvol[2] = { 0.0 };
        double inplace_surfvol[2] = { 0.0 };
        double polymass = computePolymerMass(porevol, state.saturation(), state.getCellData( state.CONCENTRATION ), poly_props_.deadPoreVol());
        double polymass_adsorbed = computePolymerAdsorbed(grid_, props_, poly_props_, state, rock_comp_props_);
        double init_polymass = polymass + polymass_adsorbed;
        double tot_injected[2] = { 0.0 };
        double tot_produced[2] = { 0.0 };
        double tot_polyinj = 0.0;
        double tot_polyprod = 0.0;
        Opm::computeSaturatedVol(porevol, state.surfacevol(), init_surfvol);
        Opm::Watercut watercut;
        watercut.push(0.0, 0.0, 0.0);
        Opm::WellReport wellreport;
        std::vector<double> fractional_flows;
        std::vector<double> well_resflows_phase;
        if (wells_) {
            well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0);
            wellreport.push(props_, *wells_, state.pressure(), state.surfacevol(),
                            state.saturation(), 0.0, well_state.bhp(), well_state.perfRates());
        }
        // Report timestep and (optionally) write state to disk.
        timer.report(std::cout);
        if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
            if (output_vtk_) {
                outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
            }
            outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
        }

        initial_pressure = state.pressure();

        // Solve pressure equation.
        if (check_well_controls_) {
            computeFractionalFlow(props_, poly_props_, allcells_,
                                  state.pressure(), state.temperature(), state.surfacevol(), state.saturation(),
                                  state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX ) ,
                                  fractional_flows);
            wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
        }
        bool well_control_passed = !check_well_controls_;
        int well_control_iteration = 0;
        do {
            // Run solver
            pressure_timer.start();
            psolver_.solve(timer.currentStepLength(), state, well_state);

            // Renormalize pressure if both fluids and rock are
            // incompressible, and there are no pressure
            // conditions (bcs or wells).  It is deemed sufficient
            // for now to renormalize using geometric volume
            // instead of pore volume.
            if (psolver_.singularPressure()) {
                // Compute average pressures of previous and last
                // step, and total volume.
                double av_prev_press = 0.0;
                double av_press = 0.0;
                double tot_vol = 0.0;
                const int num_cells = grid_.number_of_cells;
                for (int cell = 0; cell < num_cells; ++cell) {
                    av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell];
                    av_press      += state.pressure()[cell]*grid_.cell_volumes[cell];
                    tot_vol       += grid_.cell_volumes[cell];
                }
                // Renormalization constant
                const double ren_const = (av_prev_press - av_press)/tot_vol;
                for (int cell = 0; cell < num_cells; ++cell) {
                    state.pressure()[cell] += ren_const;
                }
                const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells;
                for (int well = 0; well < num_wells; ++well) {
                    well_state.bhp()[well] += ren_const;
                }
            }

            // Stop timer and report
//.........这里部分代码省略.........

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


    std::cout << "\n\n================    Starting main simulation loop     ===============\n"
              << std::flush;

    SimulatorReport rep;
    if (!use_deck) {
        // Simple simulation without a deck.
        PolymerInflowBasic polymer_inflow(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
                                          param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
                                          param.getDefault("poly_amount", poly_props.cMax()));
        WellsManager wells;
        SimulatorCompressiblePolymer simulator(param,
                                               *grid->c_grid(),
                                               *props,
                                               poly_props,
                                               rock_comp->isActive() ? rock_comp.get() : 0,
                                               wells,
                                               polymer_inflow,
                                               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;
        Opm::TimeMapPtr timeMap(new Opm::TimeMap(deck));
        SimulatorTimer simtimer;
        simtimer.init(timeMap);
        // Check for WPOLYMER presence in last report step to decide
        // polymer injection control type.
        const bool use_wpolymer = deck->hasKeyword("WPOLYMER");
        if (use_wpolymer) {
            if (param.has("poly_start_days")) {
                OPM_MESSAGE("Warning: Using WPOLYMER to control injection since it was found in deck. "
                        "You seem to be trying to control it via parameter poly_start_days (etc.) as well.");
            }
        }
        for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) {
            simtimer.setCurrentStepNum(reportStepIdx);

            // Report on start of report step.
            std::cout << "\n\n--------------    Starting report step " << reportStepIdx << "    --------------"
                      << "\n                  (number of remaining steps: "
                      << simtimer.numSteps() - step << ")\n\n" << std::flush;

            // Create new wells, polymer inflow controls.
            WellsManager wells(eclipseState , reportStepIdx , *grid->c_grid(), props->permeability());
            boost::scoped_ptr<PolymerInflowInterface> polymer_inflow;
            if (use_wpolymer) {
                if (wells.c_wells() == 0) {
                    OPM_THROW(std::runtime_error, "Cannot control polymer injection via WPOLYMER without wells.");
                }
                polymer_inflow.reset(new PolymerInflowFromDeck(eclipseState, *wells.c_wells(), props->numCells(), simtimer.currentStepNum()));
            } else {
                polymer_inflow.reset(new PolymerInflowBasic(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
                                                            param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
                                                            param.getDefault("poly_amount", poly_props.cMax())));
            }

            // @@@ HACK: we should really make a new well state and
            // properly transfer old well state to it every report step,
            // since number of wells may change etc.
            if (reportStepIdx == 0) {
                well_state.init(wells.c_wells(), *state);
            }

            // Create and run simulator.
            SimulatorCompressiblePolymer simulator(param,
                                                   *grid->c_grid(),
                                                   *props,
                                                   poly_props,
                                                   rock_comp->isActive() ? rock_comp.get() : 0,
                                                   wells,
                                                   *polymer_inflow,
                                                   linsolver,
                                                   grav);
            if (reportStepIdx == 0) {
                warnIfUnusedParams(param);
            }
            SimulatorReport epoch_rep = simulator.run(simtimer, *state, well_state);

            // 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);
}
catch (const std::exception &e) {
    std::cerr << "Program threw an exception: " << e.what() << "\n";
    throw;
}

    SimulatorReport SimulatorCompressibleAd::Impl::run(SimulatorTimer& timer,
                                                       BlackoilState& state,
                                                       WellState& well_state)
    {
        std::vector<double> transport_src;

        // Initialisation.
        std::vector<double> porevol;
        if (rock_comp_props_ && rock_comp_props_->isActive()) {
            computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
        } else {
            computePorevolume(grid_, props_.porosity(), porevol);
        }
        const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
        std::vector<double> initial_porevol = porevol;

        // Main simulation loop.
        Opm::time::StopWatch pressure_timer;
        double ptime = 0.0;
        Opm::time::StopWatch transport_timer;
        double ttime = 0.0;
        Opm::time::StopWatch step_timer;
        Opm::time::StopWatch total_timer;
        total_timer.start();
        double init_surfvol[2] = { 0.0 };
        double inplace_surfvol[2] = { 0.0 };
        double tot_injected[2] = { 0.0 };
        double tot_produced[2] = { 0.0 };
        Opm::computeSaturatedVol(porevol, state.surfacevol(), init_surfvol);
        Opm::Watercut watercut;
        watercut.push(0.0, 0.0, 0.0);
        Opm::WellReport wellreport;
        std::vector<double> fractional_flows;
        std::vector<double> well_resflows_phase;
        if (wells_) {
            well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0);
            wellreport.push(props_, *wells_,
                            state.pressure(), state.surfacevol(), state.saturation(),
                            0.0, well_state.bhp(), well_state.perfRates());
        }
        std::fstream tstep_os;
        if (output_) {
            std::string filename = output_dir_ + "/step_timing.param";
            tstep_os.open(filename.c_str(), std::fstream::out | std::fstream::app);
        }
        for (; !timer.done(); ++timer) {
            // Report timestep and (optionally) write state to disk.
            step_timer.start();
            timer.report(std::cout);
            if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
                if (output_vtk_) {
                    outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
                }
                outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
            }

            SimulatorReport sreport;

            // Solve pressure equation.
            if (check_well_controls_) {
                computeFractionalFlow(props_, allcells_,
                                      state.pressure(), state.surfacevol(), state.saturation(),
                                      fractional_flows);
                wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
            }
            bool well_control_passed = !check_well_controls_;
            int well_control_iteration = 0;
            do {
                // Run solver.
                pressure_timer.start();
                std::vector<double> initial_pressure = state.pressure();
                psolver_.solve(timer.currentStepLength(), state, well_state);

#if 0
                // Renormalize pressure if both fluids and rock are
                // incompressible, and there are no pressure
                // conditions (bcs or wells).  It is deemed sufficient
                // for now to renormalize using geometric volume
                // instead of pore volume.
                if (psolver_.singularPressure()) {
                    // Compute average pressures of previous and last
                    // step, and total volume.
                    double av_prev_press = 0.0;
                    double av_press = 0.0;
                    double tot_vol = 0.0;
                    const int num_cells = grid_.number_of_cells;
                    for (int cell = 0; cell < num_cells; ++cell) {
                        av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell];
                        av_press      += state.pressure()[cell]*grid_.cell_volumes[cell];
                        tot_vol       += grid_.cell_volumes[cell];
                    }
                    // Renormalization constant
                    const double ren_const = (av_prev_press - av_press)/tot_vol;
                    for (int cell = 0; cell < num_cells; ++cell) {
                        state.pressure()[cell] += ren_const;
                    }
                    const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells;
                    for (int well = 0; well < num_wells; ++well) {
                        well_state.bhp()[well] += ren_const;
                    }
//.........这里部分代码省略.........

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

}

//.........这里部分代码省略.........
                                               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();
        // Check for WPOLYMER presence in last epoch to decide
        // polymer injection control type.
        const bool use_wpolymer = deck->hasField("WPOLYMER");
        if (use_wpolymer) {
            if (param.has("poly_start_days")) {
                OPM_MESSAGE("Warning: Using WPOLYMER to control injection since it was found in deck. "
                        "You seem to be trying to control it via parameter poly_start_days (etc.) as well.");
            }
        }
        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) {
                    OPM_THROW(std::runtime_error, "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, polymer inflow controls.
            WellsManager wells(eclipseState , epoch , *grid->c_grid(), props->permeability());
            boost::scoped_ptr<PolymerInflowInterface> polymer_inflow;
            if (use_wpolymer) {
                if (wells.c_wells() == 0) {
                    OPM_THROW(std::runtime_error, "Cannot control polymer injection via WPOLYMER without wells.");
                }
                polymer_inflow.reset(new PolymerInflowFromDeck(*deck, *wells.c_wells(), props->numCells()));
            } else {
                polymer_inflow.reset(new PolymerInflowBasic(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
                                                            param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
                                                            param.getDefault("poly_amount", poly_props.cMax())));
            }

            // @@@ 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.
            SimulatorCompressiblePolymer simulator(param,
                                                   *grid->c_grid(),
                                                   *props,
                                                   poly_props,
                                                   rock_comp->isActive() ? rock_comp.get() : 0,
                                                   wells,
                                                   *polymer_inflow,
                                                   src,
                                                   bcs.c_bcs(),
                                                   linsolver,
                                                   grav);
            if (epoch == 0) {
                warnIfUnusedParams(param);
            }
            SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state);

            // 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);
}
catch (const std::exception &e) {
    std::cerr << "Program threw an exception: " << e.what() << "\n";
    throw;
}

    SimulatorReport SimulatorPolymer::Impl::run(SimulatorTimer& timer,
                                                PolymerState& state,
                                                WellState& well_state)
    {
        std::vector<double> transport_src;

        // Initialisation.
        std::vector<double> porevol;
        if (rock_comp_props_ && rock_comp_props_->isActive()) {
            computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
        } else {
            computePorevolume(grid_, props_.porosity(), porevol);
        }
        const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);


        // Main simulation loop.
        Opm::time::StopWatch pressure_timer;
        double ptime = 0.0;
        Opm::time::StopWatch transport_timer;
        double ttime = 0.0;
        Opm::time::StopWatch total_timer;
        total_timer.start();
        double init_satvol[2] = { 0.0 };
        double init_polymass = 0.0;
        double satvol[2] = { 0.0 };
        double polymass = 0.0;
        double polymass_adsorbed = 0.0;
        double injected[2] = { 0.0 };
        double produced[2] = { 0.0 };
        double polyinj = 0.0;
        double polyprod = 0.0;
        double tot_injected[2] = { 0.0 };
        double tot_produced[2] = { 0.0 };
        double tot_polyinj = 0.0;
        double tot_polyprod = 0.0;
        Opm::computeSaturatedVol(porevol, state.saturation(), init_satvol);
        std::cout << "\nInitial saturations are    " << init_satvol[0]/tot_porevol_init
                  << "    " << init_satvol[1]/tot_porevol_init << std::endl;
        Opm::Watercut watercut;
        watercut.push(0.0, 0.0, 0.0);
        Opm::WellReport wellreport;
        std::vector<double> fractional_flows;
        std::vector<double> well_resflows_phase;
        if (wells_) {
            well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0);
            wellreport.push(props_, *wells_, state.saturation(), 0.0, well_state.bhp(), well_state.perfRates());
        }
        for (; !timer.done(); ++timer) {
            // Report timestep and (optionally) write state to disk.
            timer.report(std::cout);
            if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
                outputState(grid_, state, timer.currentStepNum(), output_dir_);
            }

            // Solve pressure.
            do {
                pressure_timer.start();
                psolver_.solve(timer.currentStepLength(), state, well_state);
                pressure_timer.stop();
                double pt = pressure_timer.secsSinceStart();
                std::cout << "Pressure solver took:  " << pt << " seconds." << std::endl;
                ptime += pt;
            } while (false);

            // Update pore volumes if rock is compressible.
            if (rock_comp_props_ && rock_comp_props_->isActive()) {
                computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
            }

            // Process transport sources (to include bdy terms and well flows).
            Opm::computeTransportSource(grid_, src_, state.faceflux(), 1.0,
                                        wells_, well_state.perfRates(), transport_src);

            // Find inflow rate.
            const double current_time = timer.currentTime();
            double stepsize = timer.currentStepLength();
            const double inflowc0 = poly_inflow_(current_time + 1e-5*stepsize);
            const double inflowc1 = poly_inflow_(current_time + (1.0 - 1e-5)*stepsize);
            if (inflowc0 != inflowc1) {
                std::cout << "**** Warning: polymer inflow rate changes during timestep. Using rate near start of step.";
            }
            const double inflow_c = inflowc0;

            // Solve transport.
            transport_timer.start();
            if (num_transport_substeps_ != 1) {
                stepsize /= double(num_transport_substeps_);
                std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl;
            }
            for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
                tsolver_.solve(&state.faceflux()[0], &porevol[0], &transport_src[0], stepsize, inflow_c,
                               state.saturation(), state.concentration(), state.maxconcentration());
                Opm::computeInjectedProduced(props_, poly_props_,
                                             state.saturation(), state.concentration(), state.maxconcentration(),
                                             transport_src, timer.currentStepLength(), inflow_c,
                                             injected, produced, polyinj, polyprod);
                if (use_segregation_split_) {
                    tsolver_.solveGravity(columns_, &porevol[0], stepsize,
                                          state.saturation(), state.concentration(), state.maxconcentration());
//.........这里部分代码省略.........

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

    SimulatorReport rep;
    if (!use_deck) {
        std::cout << "\n\n================    Starting main simulation loop     ===============\n"
                  << "                        (number of report steps: 1)\n\n" << std::flush;
        // Simple simulation without a deck.
        WellsManager wells; // no wells.
        SimulatorIncompTwophase 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.
        Opm::TimeMapConstPtr timeMap = eclipseState->getSchedule()->getTimeMap();

        std::cout << "\n\n================    Starting main simulation loop     ===============\n"
                  << "                        (number of report steps: "
                  << timeMap->numTimesteps() << ")\n\n" << std::flush;
        WellState well_state;
        int step = 0;
        SimulatorTimer simtimer;
        // Use timer for last epoch to obtain total time.
        simtimer.init(timeMap);
        const double total_time = simtimer.totalTime();
        // for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) {
        size_t reportStepIdx = 0; // Only handle a single, unchanging well setup.
        {
            // Update the timer.
            simtimer.setCurrentStepNum(step);
            simtimer.setTotalTime(total_time);

            // Report on start of report step.
            // std::cout << "\n\n--------------    Starting report step " << reportStepIdx << "    --------------"
            //           << "\n                  (number of time steps: "
            //           << simtimer.numSteps() - step << ")\n\n" << std::flush;

            // Create new wells, well_state
            WellsManager wells(eclipseState , reportStepIdx , *grid->c_grid(), props->permeability());
            // @@@ HACK: we should really make a new well state and
            // properly transfer old well state to it every report step,
            // since number of wells may change etc.
            if (reportStepIdx == 0) {
                well_state.init(wells.c_wells(), *state);
            }

            // Create and run simulator.
            SimulatorIncompTwophase simulator(param,
                                              *grid->c_grid(),
                                              *props,
                                              rock_comp->isActive() ? rock_comp.get() : 0,
                                              wells,
                                              src,
                                              bcs.c_bcs(),
                                              linsolver,
                                              grav);
            if (reportStepIdx == 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);
    }

}
catch (const std::exception &e) {
    std::cerr << "Program threw an exception: " << e.what() << "\n";
    throw;
}

    SimulatorReport SimulatorFullyImplicitCompressiblePolymer::Impl::run(SimulatorTimer& timer,
                                                                         PolymerBlackoilState& state)
    {
        WellStateFullyImplicitBlackoil prev_well_state;
        // Initialisation.
        std::vector<double> porevol;
        if (rock_comp_props_ && rock_comp_props_->isActive()) {
            computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
        } else {
            computePorevolume(grid_, props_.porosity(), porevol);
        }
        std::vector<double> initial_porevol = porevol;

        std::vector<double> polymer_inflow_c(grid_.number_of_cells);
        // Main simulation loop.
        Opm::time::StopWatch solver_timer;
        double stime = 0.0;
        Opm::time::StopWatch step_timer;
        Opm::time::StopWatch total_timer;
        total_timer.start();
        std::string tstep_filename = output_dir_ + "/step_timing.txt";
        std::ofstream tstep_os(tstep_filename.c_str());

        //Main simulation loop.
        while (!timer.done()) {
#if 0
            double tot_injected[2] = { 0.0 };
            double tot_produced[2] = { 0.0 };
            Opm::Watercut watercut;
            watercut.push(0.0, 0.0, 0.0);
            std::vector<double> fractional_flows;
            std::vector<double> well_resflows_phase;
            if (wells_) {
                well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0);
            }
            std::fstream tstep_os;
            if (output_) {
                std::string filename = output_dir_ + "/step_timing.param";
                tstep_os.open(filename.c_str(), std::fstream::out | std::fstream::app);
            }
#endif
            // Report timestep and (optionally) write state to disk.

            step_timer.start();
            timer.report(std::cout);

            WellsManager wells_manager(eclipse_state_,
                                       timer.currentStepNum(),
                                       Opm::UgGridHelpers::numCells(grid_),
                                       Opm::UgGridHelpers::globalCell(grid_),
                                       Opm::UgGridHelpers::cartDims(grid_),
                                       Opm::UgGridHelpers::dimensions(grid_),
                                       Opm::UgGridHelpers::cell2Faces(grid_),
                                       Opm::UgGridHelpers::beginFaceCentroids(grid_),
                                       props_.permeability());
            const Wells* wells = wells_manager.c_wells();
            WellStateFullyImplicitBlackoil well_state;
            well_state.init(wells, state.blackoilState(), prev_well_state);
            //Compute polymer inflow.
            std::unique_ptr<PolymerInflowInterface> polymer_inflow_ptr;
            if (deck_->hasKeyword("WPOLYMER")) {
                if (wells_manager.c_wells() == 0) {
                    OPM_THROW(std::runtime_error, "Cannot control polymer injection via WPOLYMER without wells.");
                }
                polymer_inflow_ptr.reset(new PolymerInflowFromDeck(deck_, eclipse_state_, *wells, Opm::UgGridHelpers::numCells(grid_), timer.currentStepNum()));
            } else {
                polymer_inflow_ptr.reset(new PolymerInflowBasic(0.0*Opm::unit::day,
                                                                1.0*Opm::unit::day,
                                                                0.0));
            }
            std::vector<double> polymer_inflow_c(Opm::UgGridHelpers::numCells(grid_));
            polymer_inflow_ptr->getInflowValues(timer.simulationTimeElapsed(),
                                                timer.simulationTimeElapsed() + timer.currentStepLength(),
                                                polymer_inflow_c);

            if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
                if (output_vtk_) {
                    outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
                }
                outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
            }
            if (output_) {
                if (timer.currentStepNum() == 0) {
                    output_writer_.writeInit(timer);
                }
                output_writer_.writeTimeStep(timer, state.blackoilState(), well_state);
            }
            // Run solver.
            solver_timer.start();
            FullyImplicitCompressiblePolymerSolver solver(grid_, props_, geo_, rock_comp_props_, polymer_props_, *wells_manager.c_wells(), linsolver_);
            solver.step(timer.currentStepLength(), state, well_state, polymer_inflow_c);
            // Stop timer and report.
            solver_timer.stop();
            const double st = solver_timer.secsSinceStart();
            std::cout << "Fully implicit solver took:  " << st << " seconds." << std::endl;

            stime += st;
            // Update pore volumes if rock is compressible.
            if (rock_comp_props_ && rock_comp_props_->isActive()) {
                initial_porevol = porevol;
//.........这里部分代码省略.........

    SimulatorReport SimulatorFullyImplicitBlackoil::Impl::run(SimulatorTimer& timer,
                                                              BlackoilState& state,
                                                              WellState& well_state)
    {
        // Initialisation.
        std::vector<double> porevol;
        if (rock_comp_props_ && rock_comp_props_->isActive()) {
            computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
        } else {
            computePorevolume(grid_, props_.porosity(), porevol);
        }
        // const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
        std::vector<double> initial_porevol = porevol;

        // Main simulation loop.
        Opm::time::StopWatch solver_timer;
        double stime = 0.0;
        Opm::time::StopWatch step_timer;
        Opm::time::StopWatch total_timer;
        total_timer.start();
#if 0
        // These must be changed for three-phase.
        double init_surfvol[2] = { 0.0 };
        double inplace_surfvol[2] = { 0.0 };
        double tot_injected[2] = { 0.0 };
        double tot_produced[2] = { 0.0 };
        Opm::computeSaturatedVol(porevol, state.surfacevol(), init_surfvol);
        Opm::Watercut watercut;
        watercut.push(0.0, 0.0, 0.0);
        Opm::WellReport wellreport;
#endif
        std::vector<double> fractional_flows;
        std::vector<double> well_resflows_phase;
        if (wells_) {
            well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0);
#if 0
            wellreport.push(props_, *wells_,
                            state.pressure(), state.surfacevol(), state.saturation(),
                            0.0, well_state.bhp(), well_state.perfRates());
#endif
        }
        std::fstream tstep_os;
        if (output_) {
            std::string filename = output_dir_ + "/step_timing.param";
            tstep_os.open(filename.c_str(), std::fstream::out | std::fstream::app);
        }
        for (; !timer.done(); ++timer) {
            // Report timestep and (optionally) write state to disk.
            step_timer.start();
            timer.report(std::cout);
            if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
                if (output_vtk_) {
                    outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
                }
                outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
                outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_);

            }

            SimulatorReport sreport;

            // Solve pressure equation.
            // if (check_well_controls_) {
            //     computeFractionalFlow(props_, allcells_,
            //                           state.pressure(), state.surfacevol(), state.saturation(),
            //                           fractional_flows);
            //     wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
            // }
            bool well_control_passed = !check_well_controls_;
            int well_control_iteration = 0;
            do {
                // Run solver.
                solver_timer.start();
                std::vector<double> initial_pressure = state.pressure();
                solver_.step(timer.currentStepLength(), state, well_state);

                // Stop timer and report.
                solver_timer.stop();
                const double st = solver_timer.secsSinceStart();
                std::cout << "Fully implicit solver took:  " << st << " seconds." << std::endl;
                stime += st;
                sreport.pressure_time = st;

                // Optionally, check if well controls are satisfied.
                if (check_well_controls_) {
                    Opm::computePhaseFlowRatesPerWell(*wells_,
                                                      well_state.perfRates(),
                                                      fractional_flows,
                                                      well_resflows_phase);
                    std::cout << "Checking well conditions." << std::endl;
                    // For testing we set surface := reservoir
                    well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase);
                    ++well_control_iteration;
                    if (!well_control_passed && well_control_iteration > max_well_control_iterations_) {
                        OPM_THROW(std::runtime_error, "Could not satisfy well conditions in " << max_well_control_iterations_ << " tries.");
                    }
                    if (!well_control_passed) {
                        std::cout << "Well controls not passed, solving again." << std::endl;
                    } else {
                        std::cout << "Well conditions met." << std::endl;
//.........这里部分代码省略.........

    SimulatorReport SimulatorBase<Implementation>::run(SimulatorTimer& timer,
                                                       ReservoirState& state)
    {
        WellState prev_well_state;

        // Create timers and file for writing timing info.
        Opm::time::StopWatch solver_timer;
        double stime = 0.0;
        Opm::time::StopWatch step_timer;
        Opm::time::StopWatch total_timer;
        total_timer.start();
        std::string tstep_filename = output_writer_.outputDirectory() + "/step_timing.txt";
        std::ofstream tstep_os(tstep_filename.c_str());

        // adaptive time stepping
        std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
        if( param_.getDefault("timestep.adaptive", true ) )
        {
            adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, solver_.parallelInformation() ) );
        }

        // init output writer
        output_writer_.writeInit( timer );

        std::string restorefilename = param_.getDefault("restorefile", std::string("") );
        if( ! restorefilename.empty() )
        {
            // -1 means that we'll take the last report step that was written
            const int desiredRestoreStep = param_.getDefault("restorestep", int(-1) );
            output_writer_.restore( timer, state, prev_well_state, restorefilename, desiredRestoreStep );
        }

        unsigned int totalNewtonIterations = 0;
        unsigned int totalLinearIterations = 0;

        // Main simulation loop.
        while (!timer.done()) {
            // Report timestep.
            step_timer.start();
            if ( terminal_output_ )
            {
                timer.report(std::cout);
            }

            // Create wells and well state.
            WellsManager wells_manager(eclipse_state_,
                                       timer.currentStepNum(),
                                       Opm::UgGridHelpers::numCells(grid_),
                                       Opm::UgGridHelpers::globalCell(grid_),
                                       Opm::UgGridHelpers::cartDims(grid_),
                                       Opm::UgGridHelpers::dimensions(grid_),
                                       Opm::UgGridHelpers::cell2Faces(grid_),
                                       Opm::UgGridHelpers::beginFaceCentroids(grid_),
                                       props_.permeability(),
                                       is_parallel_run_);
            const Wells* wells = wells_manager.c_wells();
            WellState well_state;
            well_state.init(wells, state, prev_well_state);

            // give the polymer and surfactant simulators the chance to do their stuff
            asImpl().handleAdditionalWellInflow(timer, wells_manager, well_state, wells);

            // write simulation state at the report stage
            output_writer_.writeTimeStep( timer, state, well_state );

            // Max oil saturation (for VPPARS), hysteresis update.
            props_.updateSatOilMax(state.saturation());
            props_.updateSatHyst(state.saturation(), allcells_);

            // Compute reservoir volumes for RESV controls.
            asImpl().computeRESV(timer.currentStepNum(), wells, state, well_state);

            // Run a multiple steps of the solver depending on the time step control.
            solver_timer.start();

            auto solver = asImpl().createSolver(wells);

            // If sub stepping is enabled allow the solver to sub cycle
            // in case the report steps are too large for the solver to converge
            //
            // \Note: The report steps are met in any case
            // \Note: The sub stepping will require a copy of the state variables
            if( adaptiveTimeStepping ) {
                adaptiveTimeStepping->step( timer, *solver, state, well_state,  output_writer_ );
            }
            else {
                // solve for complete report step
                solver->step(timer.currentStepLength(), state, well_state);
            }

            // take time that was used to solve system for this reportStep
            solver_timer.stop();

            // accumulate the number of Newton and Linear Iterations
            totalNewtonIterations += solver->newtonIterations();
            totalLinearIterations += solver->linearIterations();

            // Report timing.
            const double st = solver_timer.secsSinceStart();

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

//.........这里部分代码省略.........
    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 (...) {
            OPM_THROW(std::runtime_error, "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: "
              << (deck->numberOfEpochs()) << ")\n\n" << std::flush;

    SimulatorReport rep;
    // 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) {
                OPM_THROW(std::runtime_error, "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.
        SimulatorFullyImplicitBlackoil simulator(param,
                                                 *grid->c_grid(),
                                                 *new_props,
                                                 rock_comp->isActive() ? rock_comp.get() : 0,
                                                 wells,
                                                 linsolver,
                                                 grav,
                                                 outputWriter);
        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);
    }

}
catch (const std::exception &e) {
    std::cerr << "Program threw an exception: " << e.what() << "\n";
    throw;
}