C++ TwophaseState::faceflux方法代码示例

    /// Compute the output.
    void IncompTpfa::computeResults(TwophaseState& state,
                                    WellState& well_state) const
    {
        // Make sure h_ contains the direct-solution matrix
        // and right hand side (not jacobian and residual).
        // TODO: optimize by only adjusting b and diagonal of A.
        UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
        ifs_tpfa_assemble(gg, &forces_, &trans_[0], &gpress_omegaweighted_[0], h_);


        // Make sure h_->x contains the direct solution vector.
        ASSERT(int(state.pressure().size()) == grid_.number_of_cells);
        ASSERT(int(state.faceflux().size()) == grid_.number_of_faces);
        std::copy(state.pressure().begin(), state.pressure().end(), h_->x);
        std::copy(well_state.bhp().begin(), well_state.bhp().end(), h_->x + grid_.number_of_cells);

        // Obtain solution.
        ifs_tpfa_solution soln = { NULL, NULL, NULL, NULL };
        soln.cell_press = &state.pressure()[0];
        soln.face_flux  = &state.faceflux()[0];
        if (wells_ != NULL) {
            ASSERT(int(well_state.bhp().size()) == wells_->number_of_wells);
            ASSERT(int(well_state.perfRates().size()) == wells_->well_connpos[ wells_->number_of_wells ]);
            soln.well_flux = &well_state.perfRates()[0];
            soln.well_press = &well_state.bhp()[0];
        }
        ifs_tpfa_press_flux(gg, &forces_, &trans_[0], h_, &soln); // TODO: Check what parts of h_ are used here.
    }

    // Solve with no rock compressibility (linear eqn).
    void IncompTpfa::solveIncomp(const double dt,
                                 TwophaseState& state,
                                 WellState& well_state)
    {
        // Set up properties.
        computePerSolveDynamicData(dt, state, well_state);

        // Assemble.
        UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
        int ok = ifs_tpfa_assemble(gg, &forces_, &trans_[0], &gpress_omegaweighted_[0], h_);
        if (!ok) {
            THROW("Failed assembling pressure system.");
        }

        // Solve.
        linsolver_.solve(h_->A, h_->b, h_->x);

        // Obtain solution.
        ASSERT(int(state.pressure().size()) == grid_.number_of_cells);
        ASSERT(int(state.faceflux().size()) == grid_.number_of_faces);
        ifs_tpfa_solution soln = { NULL, NULL, NULL, NULL };
        soln.cell_press = &state.pressure()[0];
        soln.face_flux  = &state.faceflux()[0];
        if (wells_ != NULL) {
            ASSERT(int(well_state.bhp().size()) == wells_->number_of_wells);
            ASSERT(int(well_state.perfRates().size()) == wells_->well_connpos[ wells_->number_of_wells ]);
            soln.well_flux = &well_state.perfRates()[0];
            soln.well_press = &well_state.bhp()[0];
        }
        ifs_tpfa_press_flux(gg, &forces_, &trans_[0], h_, &soln);
    }

    void TransportSolverTwophaseReorder::solve(const double* porevolume,
                                               const double* source,
                                               const double dt,
                                               TwophaseState& state)
    {
        darcyflux_ = &state.faceflux()[0];
        porevolume_ = porevolume;
        source_ = source;
        dt_ = dt;
        toWaterSat(state.saturation(), saturation_);

#ifdef EXPERIMENT_GAUSS_SEIDEL
        std::vector<int> seq(grid_.number_of_cells);
        std::vector<int> comp(grid_.number_of_cells + 1);
        int ncomp;
        compute_sequence_graph(&grid_, darcyflux_,
                               &seq[0], &comp[0], &ncomp,
                               &ia_upw_[0], &ja_upw_[0]);
        const int nf = grid_.number_of_faces;
        std::vector<double> neg_darcyflux(nf);
        std::transform(darcyflux_, darcyflux_ + nf, neg_darcyflux.begin(), std::negate<double>());
        compute_sequence_graph(&grid_, &neg_darcyflux[0],
                               &seq[0], &comp[0], &ncomp,
                               &ia_downw_[0], &ja_downw_[0]);
#endif
        std::fill(reorder_iterations_.begin(),reorder_iterations_.end(),0);
        reorderAndTransport(grid_, darcyflux_);
        toBothSat(saturation_, state.saturation());
    }

//.........这里部分代码省略.........
        int pside = param.get<int>("pside");
        double pside_pressure = param.get<double>("pside_pressure");
        bcs.pressureSide(*grid->c_grid(), FlowBCManager::Side(pside), pside_pressure);
    }

    // Linear solver.
    LinearSolverFactory linsolver(param);

    // Pressure solver.
    Opm::IncompTpfa psolver(*grid->c_grid(), *props, 0, linsolver,
                            0.0, 0.0, 0,
                            grav, wells->c_wells(), src, bcs.c_bcs());

    // Choice of tof solver.
    bool use_dg = param.getDefault("use_dg", false);
    int dg_degree = -1;
    if (use_dg) {
        dg_degree = param.getDefault("dg_degree", 0);
    }

    // Write parameters used for later reference.
    bool output = param.getDefault("output", true);
    std::ofstream epoch_os;
    std::string output_dir;
    if (output) {
        output_dir =
            param.getDefault("output_dir", std::string("output"));
        boost::filesystem::path fpath(output_dir);
        try {
            create_directories(fpath);
        }
        catch (...) {
            THROW("Creating directories failed: " << fpath);
        }
        std::string filename = output_dir + "/epoch_timing.param";
        epoch_os.open(filename.c_str(), std::fstream::trunc | std::fstream::out);
        // open file to clean it. The file is appended to in SimulatorTwophase
        filename = output_dir + "/step_timing.param";
        std::fstream step_os(filename.c_str(), std::fstream::trunc | std::fstream::out);
        step_os.close();
        param.writeParam(output_dir + "/simulation.param");
    }

    // Init wells.
    Opm::WellState well_state;
    well_state.init(wells->c_wells(), state);

    // Main solvers.
    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();
    std::cout << "\n\n================    Starting main solvers     ===============" << std::endl;

    // Solve pressure.
    pressure_timer.start();
    psolver.solve(1.0, state, well_state);
    pressure_timer.stop();
    double pt = pressure_timer.secsSinceStart();
    std::cout << "Pressure solver took:  " << pt << " seconds." << std::endl;
    ptime += pt;

    // Process transport sources (to include bdy terms and well flows).
    std::vector<double> transport_src;
    Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0,
                                wells->c_wells(), well_state.perfRates(), transport_src);

    // Solve time-of-flight.
    std::vector<double> tof;
    if (use_dg) {
        bool use_cvi = param.getDefault("use_cvi", false);
        Opm::TransportModelTracerTofDiscGal tofsolver(*grid->c_grid(), use_cvi);
        transport_timer.start();
        tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], dg_degree, tof);
        transport_timer.stop();
    } else {
        Opm::TransportModelTracerTof tofsolver(*grid->c_grid());
        transport_timer.start();
        tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof);
        transport_timer.stop();
    }
    double tt = transport_timer.secsSinceStart();
    std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
    ttime += tt;
    total_timer.stop();

    // Output.
    if (output) {
        std::string tof_filename = output_dir + "/tof.txt";
        std::ofstream tof_stream(tof_filename.c_str());
        std::copy(tof.begin(), tof.end(), std::ostream_iterator<double>(tof_stream, "\n"));
    }

    std::cout << "\n\n================    End of simulation     ===============\n"
              << "Total time taken: " << total_timer.secsSinceStart()
              << "\n  Pressure time:  " << ptime
              << "\n  Transport time: " << ttime << std::endl;
}