mirror of
https://github.com/OPM/opm-simulators.git
synced 2024-12-26 17:20:59 -06:00
Whitespace fixes.
It turns out I accidentally used tabs for a while, this commit fixes that for all touched files.
This commit is contained in:
parent
3c905845f9
commit
097542a527
@ -158,7 +158,7 @@ namespace Opm
|
|||||||
}
|
}
|
||||||
// ... then set the one corrresponding to this variable to identity.
|
// ... then set the one corrresponding to this variable to identity.
|
||||||
assert(blocksizes[index] == num_elem);
|
assert(blocksizes[index] == num_elem);
|
||||||
jac[index] = M(M::IdentityMatrix, val.size());
|
jac[index] = M(M::IdentityMatrix, val.size());
|
||||||
// if(val.size()>0)
|
// if(val.size()>0)
|
||||||
// {
|
// {
|
||||||
// // if val is empty the following will run into an assertion
|
// // if val is empty the following will run into an assertion
|
||||||
@ -343,7 +343,7 @@ namespace Opm
|
|||||||
// D D1 = val_.matrix().asDiagonal();
|
// D D1 = val_.matrix().asDiagonal();
|
||||||
// D D2 = rhs.val_.matrix().asDiagonal();
|
// D D2 = rhs.val_.matrix().asDiagonal();
|
||||||
M D1(val_.matrix().asDiagonal());
|
M D1(val_.matrix().asDiagonal());
|
||||||
M D2(rhs.val_.matrix().asDiagonal());
|
M D2(rhs.val_.matrix().asDiagonal());
|
||||||
for (int block = 0; block < num_blocks; ++block) {
|
for (int block = 0; block < num_blocks; ++block) {
|
||||||
assert(jac_[block].rows() == rhs.jac_[block].rows());
|
assert(jac_[block].rows() == rhs.jac_[block].rows());
|
||||||
assert(jac_[block].cols() == rhs.jac_[block].cols());
|
assert(jac_[block].cols() == rhs.jac_[block].cols());
|
||||||
@ -382,8 +382,8 @@ namespace Opm
|
|||||||
// D D2 = rhs.val_.matrix().asDiagonal();
|
// D D2 = rhs.val_.matrix().asDiagonal();
|
||||||
// D D3 = (1.0/(rhs.val_*rhs.val_)).matrix().asDiagonal();
|
// D D3 = (1.0/(rhs.val_*rhs.val_)).matrix().asDiagonal();
|
||||||
M D1(val_.matrix().asDiagonal());
|
M D1(val_.matrix().asDiagonal());
|
||||||
M D2(rhs.val_.matrix().asDiagonal());
|
M D2(rhs.val_.matrix().asDiagonal());
|
||||||
M D3((1.0/(rhs.val_*rhs.val_)).matrix().asDiagonal());
|
M D3((1.0/(rhs.val_*rhs.val_)).matrix().asDiagonal());
|
||||||
for (int block = 0; block < num_blocks; ++block) {
|
for (int block = 0; block < num_blocks; ++block) {
|
||||||
assert(jac_[block].rows() == rhs.jac_[block].rows());
|
assert(jac_[block].rows() == rhs.jac_[block].rows());
|
||||||
assert(jac_[block].cols() == rhs.jac_[block].cols());
|
assert(jac_[block].cols() == rhs.jac_[block].cols());
|
||||||
@ -412,8 +412,8 @@ namespace Opm
|
|||||||
int num_blocks = jac_.size();
|
int num_blocks = jac_.size();
|
||||||
os << "Value =\n" << val_ << "\n\nJacobian =\n";
|
os << "Value =\n" << val_ << "\n\nJacobian =\n";
|
||||||
for (int i = 0; i < num_blocks; ++i) {
|
for (int i = 0; i < num_blocks; ++i) {
|
||||||
Eigen::SparseMatrix<double> m;
|
Eigen::SparseMatrix<double> m;
|
||||||
jac_[i].toSparse(m);
|
jac_[i].toSparse(m);
|
||||||
os << "Sub Jacobian #" << i << '\n' << m << "\n";
|
os << "Sub Jacobian #" << i << '\n' << m << "\n";
|
||||||
}
|
}
|
||||||
return os;
|
return os;
|
||||||
|
@ -251,7 +251,7 @@ struct HelperOps
|
|||||||
}
|
}
|
||||||
|
|
||||||
private:
|
private:
|
||||||
Eigen::SparseMatrix<double> select_;
|
Eigen::SparseMatrix<double> select_;
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
||||||
@ -460,15 +460,15 @@ collapseJacs(const AutoDiffBlock<double>& x, Matrix& jacobian)
|
|||||||
int block_col_start = 0;
|
int block_col_start = 0;
|
||||||
for (int block = 0; block < nb; ++block) {
|
for (int block = 0; block < nb; ++block) {
|
||||||
const ADB::M& jac1 = x.derivative()[block];
|
const ADB::M& jac1 = x.derivative()[block];
|
||||||
Eigen::SparseMatrix<double> jac;
|
Eigen::SparseMatrix<double> jac;
|
||||||
jac1.toSparse(jac);
|
jac1.toSparse(jac);
|
||||||
for (Eigen::SparseMatrix<double>::Index k = 0; k < jac.outerSize(); ++k) {
|
for (Eigen::SparseMatrix<double>::Index k = 0; k < jac.outerSize(); ++k) {
|
||||||
for (Eigen::SparseMatrix<double>::InnerIterator i(jac, k); i ; ++i) {
|
for (Eigen::SparseMatrix<double>::InnerIterator i(jac, k); i ; ++i) {
|
||||||
if (i.value() != 0.0) {
|
if (i.value() != 0.0) {
|
||||||
t.push_back(Tri(i.row(),
|
t.push_back(Tri(i.row(),
|
||||||
i.col() + block_col_start,
|
i.col() + block_col_start,
|
||||||
i.value()));
|
i.value()));
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
block_col_start += jac.cols();
|
block_col_start += jac.cols();
|
||||||
@ -588,8 +588,8 @@ vertcatCollapseJacs(const std::vector<AutoDiffBlock<double> >& x)
|
|||||||
if (!x[elem].derivative().empty()) {
|
if (!x[elem].derivative().empty()) {
|
||||||
for (int block = 0; block < num_blocks; ++block) {
|
for (int block = 0; block < num_blocks; ++block) {
|
||||||
// const ADB::M& jac = x[elem].derivative()[block];
|
// const ADB::M& jac = x[elem].derivative()[block];
|
||||||
M jac;
|
M jac;
|
||||||
x[elem].derivative()[block].toSparse(jac);
|
x[elem].derivative()[block].toSparse(jac);
|
||||||
for (M::Index k = 0; k < jac.outerSize(); ++k) {
|
for (M::Index k = 0; k < jac.outerSize(); ++k) {
|
||||||
for (M::InnerIterator i(jac, k); i ; ++i) {
|
for (M::InnerIterator i(jac, k); i ; ++i) {
|
||||||
t.push_back(Tri(i.row() + block_row_start,
|
t.push_back(Tri(i.row() + block_row_start,
|
||||||
|
@ -230,7 +230,7 @@ namespace Opm {
|
|||||||
WellOps(const Wells* wells);
|
WellOps(const Wells* wells);
|
||||||
// M w2p; // well -> perf (scatter)
|
// M w2p; // well -> perf (scatter)
|
||||||
// M p2w; // perf -> well (gather)
|
// M p2w; // perf -> well (gather)
|
||||||
Eigen::SparseMatrix<double> w2p; // well -> perf (scatter)
|
Eigen::SparseMatrix<double> w2p; // well -> perf (scatter)
|
||||||
Eigen::SparseMatrix<double> p2w; // perf -> well (gather)
|
Eigen::SparseMatrix<double> p2w; // perf -> well (gather)
|
||||||
};
|
};
|
||||||
|
|
||||||
|
@ -797,7 +797,7 @@ namespace detail {
|
|||||||
// get reasonable initial conditions for the wells
|
// get reasonable initial conditions for the wells
|
||||||
updateWellControls(well_state);
|
updateWellControls(well_state);
|
||||||
|
|
||||||
// Create the primary variables.
|
// Create the primary variables.
|
||||||
SolutionState state = asImpl().variableState(reservoir_state, well_state);
|
SolutionState state = asImpl().variableState(reservoir_state, well_state);
|
||||||
|
|
||||||
if (initial_assembly) {
|
if (initial_assembly) {
|
||||||
@ -978,7 +978,7 @@ namespace detail {
|
|||||||
selectInjectingPerforations[c] = 1;
|
selectInjectingPerforations[c] = 1;
|
||||||
else
|
else
|
||||||
selectProducingPerforations[c] = 1;
|
selectProducingPerforations[c] = 1;
|
||||||
}
|
}
|
||||||
|
|
||||||
// HANDLE FLOW INTO WELLBORE
|
// HANDLE FLOW INTO WELLBORE
|
||||||
// compute phase volumetric rates at standard conditions
|
// compute phase volumetric rates at standard conditions
|
||||||
@ -1441,9 +1441,9 @@ namespace detail {
|
|||||||
eqs.push_back(residual_.well_eq);
|
eqs.push_back(residual_.well_eq);
|
||||||
ADB total_residual = vertcatCollapseJacs(eqs);
|
ADB total_residual = vertcatCollapseJacs(eqs);
|
||||||
const std::vector<M>& Jn = total_residual.derivative();
|
const std::vector<M>& Jn = total_residual.derivative();
|
||||||
typedef Eigen::SparseMatrix<double> Sp;
|
typedef Eigen::SparseMatrix<double> Sp;
|
||||||
Sp Jn0;
|
Sp Jn0;
|
||||||
Jn[0].toSparse(Jn0);
|
Jn[0].toSparse(Jn0);
|
||||||
const Eigen::SparseLU< Sp > solver(Jn0);
|
const Eigen::SparseLU< Sp > solver(Jn0);
|
||||||
const Eigen::VectorXd& dx = solver.solve(total_residual.value().matrix());
|
const Eigen::VectorXd& dx = solver.solve(total_residual.value().matrix());
|
||||||
assert(dx.size() == (well_state.numWells() * (well_state.numPhases()+1)));
|
assert(dx.size() == (well_state.numWells() * (well_state.numPhases()+1)));
|
||||||
@ -1520,7 +1520,7 @@ namespace detail {
|
|||||||
//Target vars
|
//Target vars
|
||||||
ADB::V bhp_targets = ADB::V::Zero(nw);
|
ADB::V bhp_targets = ADB::V::Zero(nw);
|
||||||
ADB::V rate_targets = ADB::V::Zero(nw);
|
ADB::V rate_targets = ADB::V::Zero(nw);
|
||||||
Eigen::SparseMatrix<double> rate_distr(nw, np*nw);
|
Eigen::SparseMatrix<double> rate_distr(nw, np*nw);
|
||||||
|
|
||||||
//Selection variables
|
//Selection variables
|
||||||
std::vector<int> bhp_elems;
|
std::vector<int> bhp_elems;
|
||||||
@ -1633,7 +1633,7 @@ namespace detail {
|
|||||||
// For wells that are dead (not flowing), and therefore not communicating
|
// For wells that are dead (not flowing), and therefore not communicating
|
||||||
// with the reservoir, we set the equation to be equal to the well's total
|
// with the reservoir, we set the equation to be equal to the well's total
|
||||||
// flow. This will be a solution only if the target rate is also zero.
|
// flow. This will be a solution only if the target rate is also zero.
|
||||||
Eigen::SparseMatrix<double> rate_summer(nw, np*nw);
|
Eigen::SparseMatrix<double> rate_summer(nw, np*nw);
|
||||||
for (int w = 0; w < nw; ++w) {
|
for (int w = 0; w < nw; ++w) {
|
||||||
for (int phase = 0; phase < np; ++phase) {
|
for (int phase = 0; phase < np; ++phase) {
|
||||||
rate_summer.insert(w, phase*nw + w) = 1.0;
|
rate_summer.insert(w, phase*nw + w) = 1.0;
|
||||||
|
@ -432,7 +432,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
|
|||||||
ADB::M temp;
|
ADB::M temp;
|
||||||
fastSparseProduct(dmudr_diag, rs.derivative()[block], temp);
|
fastSparseProduct(dmudr_diag, rs.derivative()[block], temp);
|
||||||
// jacs[block] += temp;
|
// jacs[block] += temp;
|
||||||
jacs[block] = jacs[block] + temp;
|
jacs[block] = jacs[block] + temp;
|
||||||
}
|
}
|
||||||
return ADB::function(std::move(mu), std::move(jacs));
|
return ADB::function(std::move(mu), std::move(jacs));
|
||||||
}
|
}
|
||||||
@ -474,7 +474,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
|
|||||||
ADB::M temp;
|
ADB::M temp;
|
||||||
fastSparseProduct(dmudr_diag, rv.derivative()[block], temp);
|
fastSparseProduct(dmudr_diag, rv.derivative()[block], temp);
|
||||||
// jacs[block] += temp;
|
// jacs[block] += temp;
|
||||||
jacs[block] = jacs[block] + temp;
|
jacs[block] = jacs[block] + temp;
|
||||||
}
|
}
|
||||||
return ADB::function(std::move(mu), std::move(jacs));
|
return ADB::function(std::move(mu), std::move(jacs));
|
||||||
}
|
}
|
||||||
@ -555,7 +555,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
|
|||||||
ADB::M temp;
|
ADB::M temp;
|
||||||
fastSparseProduct(dbdr_diag, rs.derivative()[block], temp);
|
fastSparseProduct(dbdr_diag, rs.derivative()[block], temp);
|
||||||
// jacs[block] += temp;
|
// jacs[block] += temp;
|
||||||
jacs[block] = jacs[block] + temp;
|
jacs[block] = jacs[block] + temp;
|
||||||
}
|
}
|
||||||
return ADB::function(std::move(b), std::move(jacs));
|
return ADB::function(std::move(b), std::move(jacs));
|
||||||
}
|
}
|
||||||
@ -598,7 +598,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
|
|||||||
ADB::M temp;
|
ADB::M temp;
|
||||||
fastSparseProduct(dbdr_diag, rv.derivative()[block], temp);
|
fastSparseProduct(dbdr_diag, rv.derivative()[block], temp);
|
||||||
// jacs[block] += temp;
|
// jacs[block] += temp;
|
||||||
jacs[block] = jacs[block] + temp;
|
jacs[block] = jacs[block] + temp;
|
||||||
}
|
}
|
||||||
return ADB::function(std::move(b), std::move(jacs));
|
return ADB::function(std::move(b), std::move(jacs));
|
||||||
}
|
}
|
||||||
@ -822,7 +822,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
|
|||||||
}
|
}
|
||||||
return adbCapPressures;
|
return adbCapPressures;
|
||||||
}
|
}
|
||||||
|
|
||||||
/// Saturation update for hysteresis behavior.
|
/// Saturation update for hysteresis behavior.
|
||||||
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
|
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
|
||||||
void BlackoilPropsAdFromDeck::updateSatHyst(const std::vector<double>& saturation,
|
void BlackoilPropsAdFromDeck::updateSatHyst(const std::vector<double>& saturation,
|
||||||
@ -831,7 +831,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
|
|||||||
const int n = cells.size();
|
const int n = cells.size();
|
||||||
satprops_->updateSatHyst(n, cells.data(), saturation.data());
|
satprops_->updateSatHyst(n, cells.data(), saturation.data());
|
||||||
}
|
}
|
||||||
|
|
||||||
/// Update for max oil saturation.
|
/// Update for max oil saturation.
|
||||||
void BlackoilPropsAdFromDeck::updateSatOilMax(const std::vector<double>& saturation)
|
void BlackoilPropsAdFromDeck::updateSatOilMax(const std::vector<double>& saturation)
|
||||||
{
|
{
|
||||||
@ -863,7 +863,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
/// Apply correction to rs/rv according to kw VAPPARS
|
/// Apply correction to rs/rv according to kw VAPPARS
|
||||||
/// \param[in/out] r Array of n rs/rv values.
|
/// \param[in/out] r Array of n rs/rv values.
|
||||||
/// \param[in] so Array of n oil saturation values.
|
/// \param[in] so Array of n oil saturation values.
|
||||||
@ -874,7 +874,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
|
|||||||
const std::vector<int>& cells,
|
const std::vector<int>& cells,
|
||||||
const double vap) const
|
const double vap) const
|
||||||
{
|
{
|
||||||
if (!satOilMax_.empty() && vap > 0.0) {
|
if (!satOilMax_.empty() && vap > 0.0) {
|
||||||
const int n = cells.size();
|
const int n = cells.size();
|
||||||
V factor = V::Ones(n, 1);
|
V factor = V::Ones(n, 1);
|
||||||
const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
|
const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
|
||||||
@ -883,12 +883,12 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
|
|||||||
// guard against too small saturation values.
|
// guard against too small saturation values.
|
||||||
const double so_i= std::max(so[i],eps_sqrt);
|
const double so_i= std::max(so[i],eps_sqrt);
|
||||||
factor[i] = std::pow(so_i/satOilMax_[cells[i]], vap);
|
factor[i] = std::pow(so_i/satOilMax_[cells[i]], vap);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
r = factor*r;
|
r = factor*r;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
/// Apply correction to rs/rv according to kw VAPPARS
|
/// Apply correction to rs/rv according to kw VAPPARS
|
||||||
/// \param[in/out] r Array of n rs/rv values.
|
/// \param[in/out] r Array of n rs/rv values.
|
||||||
/// \param[in] so Array of n oil saturation values.
|
/// \param[in] so Array of n oil saturation values.
|
||||||
@ -899,7 +899,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
|
|||||||
const std::vector<int>& cells,
|
const std::vector<int>& cells,
|
||||||
const double vap) const
|
const double vap) const
|
||||||
{
|
{
|
||||||
if (!satOilMax_.empty() && vap > 0.0) {
|
if (!satOilMax_.empty() && vap > 0.0) {
|
||||||
const int n = cells.size();
|
const int n = cells.size();
|
||||||
V factor = V::Ones(n, 1);
|
V factor = V::Ones(n, 1);
|
||||||
const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
|
const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
|
||||||
|
@ -208,13 +208,13 @@ namespace Opm
|
|||||||
// Find sparsity structure as union of basic block sparsity structures,
|
// Find sparsity structure as union of basic block sparsity structures,
|
||||||
// corresponding to the jacobians with respect to pressure.
|
// corresponding to the jacobians with respect to pressure.
|
||||||
// Use addition to get to the union structure.
|
// Use addition to get to the union structure.
|
||||||
typedef Eigen::SparseMatrix<double> Sp;
|
typedef Eigen::SparseMatrix<double> Sp;
|
||||||
Sp structure;
|
Sp structure;
|
||||||
eqs[0].derivative()[0].toSparse(structure);
|
eqs[0].derivative()[0].toSparse(structure);
|
||||||
for (int phase = 1; phase < np; ++phase) {
|
for (int phase = 1; phase < np; ++phase) {
|
||||||
Sp s0;
|
Sp s0;
|
||||||
eqs[phase].derivative()[0].toSparse(s0);
|
eqs[phase].derivative()[0].toSparse(s0);
|
||||||
structure += s0;
|
structure += s0;
|
||||||
}
|
}
|
||||||
|
|
||||||
Eigen::SparseMatrix<double, Eigen::RowMajor> s = structure;
|
Eigen::SparseMatrix<double, Eigen::RowMajor> s = structure;
|
||||||
|
@ -47,7 +47,7 @@ namespace Opm
|
|||||||
NewtonIterationBlackoilSimple::SolutionVector
|
NewtonIterationBlackoilSimple::SolutionVector
|
||||||
NewtonIterationBlackoilSimple::computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const
|
NewtonIterationBlackoilSimple::computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const
|
||||||
{
|
{
|
||||||
/*
|
/*
|
||||||
typedef LinearisedBlackoilResidual::ADB ADB;
|
typedef LinearisedBlackoilResidual::ADB ADB;
|
||||||
const int np = residual.material_balance_eq.size();
|
const int np = residual.material_balance_eq.size();
|
||||||
ADB mass_res = residual.material_balance_eq[0];
|
ADB mass_res = residual.material_balance_eq[0];
|
||||||
@ -74,7 +74,7 @@ namespace Opm
|
|||||||
"Linear solver convergence failure.");
|
"Linear solver convergence failure.");
|
||||||
}
|
}
|
||||||
return dx;
|
return dx;
|
||||||
*/
|
*/
|
||||||
}
|
}
|
||||||
|
|
||||||
const boost::any& NewtonIterationBlackoilSimple::parallelInformation() const
|
const boost::any& NewtonIterationBlackoilSimple::parallelInformation() const
|
||||||
|
@ -62,9 +62,9 @@ namespace Opm
|
|||||||
const std::vector<M>& Jn = eqs[n].derivative();
|
const std::vector<M>& Jn = eqs[n].derivative();
|
||||||
|
|
||||||
// Use sparse LU to solve the block submatrices i.e compute inv(D)
|
// Use sparse LU to solve the block submatrices i.e compute inv(D)
|
||||||
typedef Eigen::SparseMatrix<double> Sp;
|
typedef Eigen::SparseMatrix<double> Sp;
|
||||||
Sp Jnn;
|
Sp Jnn;
|
||||||
Jn[n].toSparse(Jnn);
|
Jn[n].toSparse(Jnn);
|
||||||
#if HAVE_UMFPACK
|
#if HAVE_UMFPACK
|
||||||
const Eigen::UmfPackLU<Sp> solver(Jnn);
|
const Eigen::UmfPackLU<Sp> solver(Jnn);
|
||||||
#else
|
#else
|
||||||
@ -109,8 +109,8 @@ namespace Opm
|
|||||||
M Bu;
|
M Bu;
|
||||||
fastSparseProduct(B, u, Bu);
|
fastSparseProduct(B, u, Bu);
|
||||||
// J -= Bu;
|
// J -= Bu;
|
||||||
Bu = Bu * -1.0;
|
Bu = Bu * -1.0;
|
||||||
J = J + Bu;
|
J = J + Bu;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -150,9 +150,9 @@ namespace Opm
|
|||||||
const M& C = eq_coll.derivative()[0];
|
const M& C = eq_coll.derivative()[0];
|
||||||
|
|
||||||
// Use sparse LU to solve the block submatrices
|
// Use sparse LU to solve the block submatrices
|
||||||
typedef Eigen::SparseMatrix<double> Sp;
|
typedef Eigen::SparseMatrix<double> Sp;
|
||||||
Sp D;
|
Sp D;
|
||||||
D1.toSparse(D);
|
D1.toSparse(D);
|
||||||
#if HAVE_UMFPACK
|
#if HAVE_UMFPACK
|
||||||
const Eigen::UmfPackLU<Sp> solver(D);
|
const Eigen::UmfPackLU<Sp> solver(D);
|
||||||
#else
|
#else
|
||||||
@ -198,7 +198,7 @@ namespace Opm
|
|||||||
Eigen::SparseMatrix<double, Eigen::RowMajor>& A,
|
Eigen::SparseMatrix<double, Eigen::RowMajor>& A,
|
||||||
V& b)
|
V& b)
|
||||||
{
|
{
|
||||||
/*
|
/*
|
||||||
if (num_phases != 3) {
|
if (num_phases != 3) {
|
||||||
OPM_THROW(std::logic_error, "formEllipticSystem() requires 3 phases.");
|
OPM_THROW(std::logic_error, "formEllipticSystem() requires 3 phases.");
|
||||||
}
|
}
|
||||||
@ -278,7 +278,7 @@ namespace Opm
|
|||||||
// Create output as product of L with equations.
|
// Create output as product of L with equations.
|
||||||
A = L * total_residual.derivative()[0];
|
A = L * total_residual.derivative()[0];
|
||||||
b = L * total_residual.value().matrix();
|
b = L * total_residual.value().matrix();
|
||||||
*/
|
*/
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
@ -144,18 +144,18 @@ void fastSparseProduct(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
|||||||
|
|
||||||
inline void fastDiagSparseProduct(// const Eigen::DiagonalMatrix<double, Eigen::Dynamic>& lhs,
|
inline void fastDiagSparseProduct(// const Eigen::DiagonalMatrix<double, Eigen::Dynamic>& lhs,
|
||||||
const std::vector<double>& lhs,
|
const std::vector<double>& lhs,
|
||||||
const Eigen::SparseMatrix<double>& rhs,
|
const Eigen::SparseMatrix<double>& rhs,
|
||||||
Eigen::SparseMatrix<double>& res)
|
Eigen::SparseMatrix<double>& res)
|
||||||
{
|
{
|
||||||
res = rhs;
|
res = rhs;
|
||||||
|
|
||||||
// Multiply rows by diagonal lhs.
|
// Multiply rows by diagonal lhs.
|
||||||
int n = res.cols();
|
int n = res.cols();
|
||||||
for (int col = 0; col < n; ++col) {
|
for (int col = 0; col < n; ++col) {
|
||||||
typedef Eigen::SparseMatrix<double>::InnerIterator It;
|
typedef Eigen::SparseMatrix<double>::InnerIterator It;
|
||||||
for (It it(res, col); it; ++it) {
|
for (It it(res, col); it; ++it) {
|
||||||
it.valueRef() *= lhs[it.row()]; // lhs.diagonal()(it.row());
|
it.valueRef() *= lhs[it.row()]; // lhs.diagonal()(it.row());
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -163,19 +163,19 @@ inline void fastDiagSparseProduct(// const Eigen::DiagonalMatrix<double, Eigen::
|
|||||||
|
|
||||||
|
|
||||||
inline void fastSparseDiagProduct(const Eigen::SparseMatrix<double>& lhs,
|
inline void fastSparseDiagProduct(const Eigen::SparseMatrix<double>& lhs,
|
||||||
// const Eigen::DiagonalMatrix<double, Eigen::Dynamic>& rhs,
|
// const Eigen::DiagonalMatrix<double, Eigen::Dynamic>& rhs,
|
||||||
const std::vector<double>& rhs,
|
const std::vector<double>& rhs,
|
||||||
Eigen::SparseMatrix<double>& res)
|
Eigen::SparseMatrix<double>& res)
|
||||||
{
|
{
|
||||||
res = lhs;
|
res = lhs;
|
||||||
|
|
||||||
// Multiply columns by diagonal rhs.
|
// Multiply columns by diagonal rhs.
|
||||||
int n = res.cols();
|
int n = res.cols();
|
||||||
for (int col = 0; col < n; ++col) {
|
for (int col = 0; col < n; ++col) {
|
||||||
typedef Eigen::SparseMatrix<double>::InnerIterator It;
|
typedef Eigen::SparseMatrix<double>::InnerIterator It;
|
||||||
for (It it(res, col); it; ++it) {
|
for (It it(res, col); it; ++it) {
|
||||||
it.valueRef() *= rhs[col]; // rhs.diagonal()(col);
|
it.valueRef() *= rhs[col]; // rhs.diagonal()(col);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -36,7 +36,7 @@ using namespace Opm;
|
|||||||
|
|
||||||
bool
|
bool
|
||||||
operator ==(const Eigen::SparseMatrix<double>& A,
|
operator ==(const Eigen::SparseMatrix<double>& A,
|
||||||
const Eigen::SparseMatrix<double>& B)
|
const Eigen::SparseMatrix<double>& B)
|
||||||
{
|
{
|
||||||
// Two SparseMatrices are equal if
|
// Two SparseMatrices are equal if
|
||||||
// 0) They have the same ordering (enforced by equal types)
|
// 0) They have the same ordering (enforced by equal types)
|
||||||
@ -53,15 +53,15 @@ operator ==(const Eigen::SparseMatrix<double>& A,
|
|||||||
eq = eq && (A.nonZeros() == B.nonZeros());
|
eq = eq && (A.nonZeros() == B.nonZeros());
|
||||||
|
|
||||||
for (typename Eigen::SparseMatrix<double>::Index
|
for (typename Eigen::SparseMatrix<double>::Index
|
||||||
k0 = 0, kend = A.outerSize(); eq && (k0 < kend); ++k0) {
|
k0 = 0, kend = A.outerSize(); eq && (k0 < kend); ++k0) {
|
||||||
for (typename Eigen::SparseMatrix<double>::InnerIterator
|
for (typename Eigen::SparseMatrix<double>::InnerIterator
|
||||||
iA(A, k0), iB(B, k0); eq && (iA && iB); ++iA, ++iB) {
|
iA(A, k0), iB(B, k0); eq && (iA && iB); ++iA, ++iB) {
|
||||||
// 3) Sparsity structure
|
// 3) Sparsity structure
|
||||||
eq = (iA.row() == iB.row()) && (iA.col() == iB.col());
|
eq = (iA.row() == iB.row()) && (iA.col() == iB.col());
|
||||||
|
|
||||||
// 4) Equal non-zero elements
|
// 4) Equal non-zero elements
|
||||||
eq = eq && (iA.value() == iB.value());
|
eq = eq && (iA.value() == iB.value());
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
return eq;
|
return eq;
|
||||||
@ -143,7 +143,7 @@ BOOST_AUTO_TEST_CASE(AdditionOps)
|
|||||||
s1 <<
|
s1 <<
|
||||||
1.0, 0.0, 2.0,
|
1.0, 0.0, 2.0,
|
||||||
0.0, 1.0, 0.0,
|
0.0, 1.0, 0.0,
|
||||||
0.0, 0.0, 2.0;
|
0.0, 0.0, 2.0;
|
||||||
Sp ss(s1.sparseView());
|
Sp ss(s1.sparseView());
|
||||||
Mat s = Mat(ss);
|
Mat s = Mat(ss);
|
||||||
|
|
||||||
@ -191,7 +191,7 @@ BOOST_AUTO_TEST_CASE(MultOps)
|
|||||||
s1 <<
|
s1 <<
|
||||||
1.0, 0.0, 2.0,
|
1.0, 0.0, 2.0,
|
||||||
0.0, 1.0, 0.0,
|
0.0, 1.0, 0.0,
|
||||||
0.0, 0.0, 2.0;
|
0.0, 0.0, 2.0;
|
||||||
Sp ss(s1.sparseView());
|
Sp ss(s1.sparseView());
|
||||||
Mat s = Mat(ss);
|
Mat s = Mat(ss);
|
||||||
|
|
||||||
|
Loading…
Reference in New Issue
Block a user