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

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

    void TransportSolverTwophaseReorder::solveGravity(const double* porevolume,
                                                      const double dt,
                                                      TwophaseState& state)
    {
        // Initialize mobilities.
        const int nc = grid_.number_of_cells;
        std::vector<int> cells(nc);
        for (int c = 0; c < nc; ++c) {
            cells[c] = c;
        }
        mob_.resize(2*nc);
        props_.relperm(cells.size(), &state.saturation()[0], &cells[0], &mob_[0], 0);
        const double* mu = props_.viscosity();
        for (int c = 0; c < nc; ++c) {
            mob_[2*c] /= mu[0];
            mob_[2*c + 1] /= mu[1];
        }

        // Set up other variables.
        porevolume_ = porevolume;
        dt_ = dt;
        toWaterSat(state.saturation(), saturation_);

        // Solve on all columns.
        int num_iters = 0;
        for (std::vector<std::vector<int> >::size_type i = 0; i < columns_.size(); i++) {
            // std::cout << "==== new column" << std::endl;
            num_iters += solveGravityColumn(columns_[i]);
        }
        std::cout << "Gauss-Seidel column solver average iterations: "
                  << double(num_iters)/double(columns_.size()) << std::endl;

        toBothSat(saturation_, state.saturation());
    }

void
VertEqImpl::downscale (const TwophaseState &coarseScale,
                       TwophaseState &fineScale) {
	// assume that the fineScale storage is already initialized
	if (!fineScale.pressure().size() == ts->number_of_cells) {
		throw OPM_EXC ("Fine scale state is not dimensioned correctly");
	}

	// properties object handle the actual downscaling since it
	// already has the information about the interface.
	// update the coarse saturation *before* we downscale to 3D,
	// since we need the residual interface for that.
	pr->upd_res_sat (&coarseScale.saturation ()[0]);
	pr->downscale_saturation (&coarseScale.saturation ()[0],
	                          &fineScale.saturation ()[0]);
	pr->downscale_pressure (&coarseScale.saturation ()[0],
	                        &coarseScale.pressure ()[0],
	                        &fineScale.pressure ()[0]);
}

    /// Compute per-solve dynamic properties.
    void IncompTpfa::computePerSolveDynamicData(const double /*dt*/,
                                                const TwophaseState& state,
                                                const WellState& /*well_state*/)
    {
        // Computed here:
        //
        // std::vector<double> wdp_;
        // std::vector<double> totmob_;
        // std::vector<double> omega_;
        // std::vector<double> trans_;
        // std::vector<double> gpress_omegaweighted_;
        // std::vector<double> initial_porevol_;
        // ifs_tpfa_forces forces_;

        // wdp_
        if (wells_) {
            Opm::computeWDP(*wells_, grid_, state.saturation(), props_.density(),
                            gravity_ ? gravity_[2] : 0.0, true, wdp_);
        }
        // totmob_, omega_, gpress_omegaweighted_
        if (gravity_) {
            computeTotalMobilityOmega(props_, allcells_, state.saturation(), totmob_, omega_);
            mim_ip_density_update(grid_.number_of_cells, grid_.cell_facepos,
                                  &omega_[0],
                                  &gpress_[0], &gpress_omegaweighted_[0]);
        } else {
            computeTotalMobility(props_, allcells_, state.saturation(), totmob_);
        }
        // trans_
        tpfa_eff_trans_compute(const_cast<UnstructuredGrid*>(&grid_), &totmob_[0], &htrans_[0], &trans_[0]);
        // initial_porevol_
        if (rock_comp_props_ && rock_comp_props_->isActive()) {
            computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), initial_porevol_);
        }
        // forces_
        forces_.src = src_.empty() ? NULL : &src_[0];
        forces_.bc = bcs_;
        forces_.W = wells_;
        forces_.totmob = &totmob_[0];
        forces_.wdp = wdp_.empty() ? NULL : &wdp_[0];
    }

void
VertEqImpl::upscale (const TwophaseState& fineScale,
                     TwophaseState& coarseScale) {
	// dimension state object to the top grid
	coarseScale.init (*ts, pr->numPhases ());

	// upscale pressure and saturation to find the initial state of
	// the two-dimensional domain. we only need to set the pressure
	// and saturation, the flux is an output field. these methods
	// are handled by the props class, since it already has access to
	// the densities and weights.
	pr->upscale_saturation (&fineScale.saturation ()[0],
	                        &coarseScale.saturation ()[0]);
	pr->upd_res_sat (&coarseScale.saturation ()[0]);
	pr->upscale_pressure (&coarseScale.saturation ()[0],
	                      &fineScale.pressure ()[0],
	                      &coarseScale.pressure ()[0]);

	// use the regular helper method to initialize the face pressure
	// since it is implemented in the header, we have access to it
	// even though it is in an anonymous namespace!
	const UnstructuredGrid& g = this->grid();

	initFacePressure (UgGridHelpers::dimensions (g),
	                  UgGridHelpers::numFaces (g),
	                  UgGridHelpers::faceCells (g),
	                  UgGridHelpers::beginFaceCentroids (g),
	                  UgGridHelpers::beginCellCentroids (g),
	                  coarseScale);

	// update the properties from the initial state (the
	// simulation object won't call this method before the
	// first timestep; it assumes that the state is initialized
	// accordingly (which is what we do here now)
	notify (coarseScale);
}

//.........这里部分代码省略.........
    /// \internal [pore volume]
    /// \endinternal

    /// \page tutorial3
    /// \details Set up the transport solver. This is a reordering implicit Euler transport solver.
    /// \snippet tutorial3.cpp transport solver
    /// \internal [transport solver]
    const double tolerance = 1e-9;
    const int max_iterations = 30;
    Opm::TransportSolverTwophaseReorder transport_solver(grid, props, NULL, tolerance, max_iterations);
    /// \internal [transport solver]
    /// \endinternal

    /// \page tutorial3
    /// \details Time integration parameters
    /// \snippet tutorial3.cpp time parameters
    /// \internal [time parameters]
    const double dt = 0.1*day;
    const int num_time_steps = 20;
    /// \internal [time parameters]
    /// \endinternal


    /// \page tutorial3
    /// \details We define a vector which contains all cell indexes. We use this
    /// vector to set up parameters on the whole domain.
    /// \snippet tutorial3.cpp cell indexes
    /// \internal [cell indexes]
    std::vector<int> allcells(num_cells);
    for (int cell = 0; cell < num_cells; ++cell) {
        allcells[cell] = cell;
    }
    /// \internal [cell indexes]
    /// \endinternal

    /// \page tutorial3
    /// \details
    /// We set up a two-phase state object, and
    /// initialize water saturation to minimum everywhere.
    /// \snippet tutorial3.cpp two-phase state
    /// \internal [two-phase state]
    TwophaseState state;
    state.init(grid.number_of_cells , grid.number_of_faces, 2);
    initSaturation( allcells , props , state , MinSat );

    /// \internal [two-phase state]
    /// \endinternal

    /// \page tutorial3
    /// \details This string stream will be used to construct a new
    /// output filename at each timestep.
    /// \snippet tutorial3.cpp output stream
    /// \internal [output stream]
    std::ostringstream vtkfilename;
    /// \internal [output stream]
    /// \endinternal


    /// \page tutorial3
    /// \details Loop over the time steps.
    /// \snippet tutorial3.cpp time loop
    /// \internal [time loop]
    for (int i = 0; i < num_time_steps; ++i) {
        /// \internal [time loop]
	/// \endinternal


        /// \page tutorial3
        /// \details Solve the pressure equation
        /// \snippet tutorial3.cpp solve pressure
        /// \internal [solve pressure]
        psolver.solve(dt, state, well_state);
        /// \internal [solve pressure]
	/// \endinternal

        /// \page tutorial3
        /// \details  Solve the transport equation.
        /// \snippet tutorial3.cpp transport solve
	/// \internal [transport solve]
        transport_solver.solve(&porevol[0], &src[0], dt, state);
        /// \internal [transport solve]
	/// \endinternal

        /// \page tutorial3
        /// \details Write the output to file.
        /// \snippet tutorial3.cpp write output
	/// \internal [write output]
        vtkfilename.str("");
        vtkfilename << "tutorial3-" << std::setw(3) << std::setfill('0') << i << ".vtu";
        std::ofstream vtkfile(vtkfilename.str().c_str());
        Opm::DataMap dm;
        dm["saturation"] = &state.saturation();
        dm["pressure"] = &state.pressure();
        Opm::writeVtkData(grid, dm, vtkfile);
    }
}
catch (const std::exception &e) {
    std::cerr << "Program threw an exception: " << e.what() << "\n";
    throw;
}