OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #233 from atgeirr/improve-cpr-performance

Browse Source 
Improve performance of CPR preconditioner.
This commit is contained in:
Bård Skaflestad 2014-11-14 08:47:05 +01:00
commit 7ac155531a
1 changed files with 21 additions and 15 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

34
opm/autodiff/NewtonIterationBlackoilCPR.cpp
View File

@ -44,7 +44,7 @@
#include <opm/core/utility/platform_dependent/reenable_warnings.h> #include <opm/core/utility/platform_dependent/reenable_warnings.h>
#include <Eigen/SparseLU> #include <Eigen/UmfPackSupport>
namespace Opm namespace Opm
{ {
@ -249,13 +249,29 @@ 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)
const Eigen::SparseLU< M > solver(Jn[n]); const Eigen::UmfPackLU< M > solver(Jn[n]);
M id(Jn[n].rows(), Jn[n].cols());
id.setIdentity();
const M Di = solver.solve(id);
// compute inv(D)*bn for the update of the right hand side // compute inv(D)*bn for the update of the right hand side
const Eigen::VectorXd& Dibn = solver.solve(eqs[n].value().matrix()); const Eigen::VectorXd& Dibn = solver.solve(eqs[n].value().matrix());
std::vector<V> vals(num_eq); // Number n will remain empty. std::vector<V> vals(num_eq); // Number n will remain empty.
std::vector<std::vector<M>> jacs(num_eq); // Number n will remain empty. std::vector<std::vector<M>> jacs(num_eq); // Number n will remain empty.
for (int eq = 0; eq < num_eq; ++eq) {
jacs[eq].reserve(num_eq - 1);
const std::vector<M>& Je = eqs[eq].derivative();
const M& B = Je[n];
// Update right hand side.
vals[eq] = eqs[eq].value().matrix() - B * Dibn;
}
for (int var = 0; var < num_eq; ++var) {
if (var == n) {
continue;
}
// solve Du = C
const M u = Di * Jn[var]; // solver.solve(Jn[var]);
for (int eq = 0; eq < num_eq; ++eq) { for (int eq = 0; eq < num_eq; ++eq) {
if (eq == n) { if (eq == n) {
continue; continue;
@ -263,23 +279,13 @@ namespace Opm
const std::vector<M>& Je = eqs[eq].derivative(); const std::vector<M>& Je = eqs[eq].derivative();
const M& B = Je[n]; const M& B = Je[n];
jacs[eq].reserve(num_eq - 1);
for (int var = 0; var < num_eq; ++var) {
if (var == n) {
continue;
}
// Create new jacobians. // Create new jacobians.
// Add A // Add A
jacs[eq].push_back(Je[var]); jacs[eq].push_back(Je[var]);
M& J = jacs[eq].back(); M& J = jacs[eq].back();
// solve Du = C // Subtract Bu (B*inv(D)*C)
const M& u = solver.solve(Jn[var]);
// Subtract Bu (B*inv(D))
J -= B * u; J -= B * u;
} }
// Update right hand side.
vals[eq] = eqs[eq].value().matrix() - B * Dibn;
} }
// Create return value. // Create return value.
@ -317,7 +323,7 @@ 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
const Eigen::SparseLU< M > solver(D); const Eigen::UmfPackLU< M > solver(D);
// Compute value of eliminated variable. // Compute value of eliminated variable.
const Eigen::VectorXd b = (equation.value().matrix() - C * partial_solution.matrix()); const Eigen::VectorXd b = (equation.value().matrix() - C * partial_solution.matrix());