mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
Merge pull request #230 from totto82/fixSchur
Solve sub matrix systems in the Schur complement
This commit is contained in:
commit
e5aef85295
@ -44,6 +44,7 @@
|
|||||||
|
|
||||||
#include <opm/core/utility/platform_dependent/reenable_warnings.h>
|
#include <opm/core/utility/platform_dependent/reenable_warnings.h>
|
||||||
|
|
||||||
|
#include <Eigen/SparseLU>
|
||||||
|
|
||||||
namespace Opm
|
namespace Opm
|
||||||
{
|
{
|
||||||
@ -193,7 +194,7 @@ namespace Opm
|
|||||||
|
|
||||||
// Construct linear solver.
|
// Construct linear solver.
|
||||||
const double tolerance = 1e-3;
|
const double tolerance = 1e-3;
|
||||||
const int maxit = 5000;
|
const int maxit = 150;
|
||||||
const int verbosity = 0;
|
const int verbosity = 0;
|
||||||
const int restart = 40;
|
const int restart = 40;
|
||||||
Dune::RestartedGMResSolver<Vector> linsolve(opA, sp, precond, tolerance, restart, maxit, verbosity);
|
Dune::RestartedGMResSolver<Vector> linsolve(opA, sp, precond, tolerance, restart, maxit, verbosity);
|
||||||
@ -240,36 +241,45 @@ namespace Opm
|
|||||||
}
|
}
|
||||||
|
|
||||||
// Schur complement of (A B ; C D) wrt. D is A - B*inv(D)*C.
|
// Schur complement of (A B ; C D) wrt. D is A - B*inv(D)*C.
|
||||||
// This is applied to all 2x2 block submatrices.
|
// This is applied to all 2x2 block submatrices
|
||||||
// We require that D is diagonal.
|
// The right hand side is modified accordingly. bi = bi - B * inv(D)* bn;
|
||||||
const M& D = eqs[n].derivative()[n];
|
// We do not explicitly compute inv(D) instead Du = C is solved
|
||||||
if (!isDiagonal(D)) {
|
|
||||||
// std::cout << "++++++++++++++++++++++++++++++++++++++++++++\n"
|
// Extract the submatrix
|
||||||
// << D
|
const std::vector<M>& Jn = eqs[n].derivative();
|
||||||
// << "++++++++++++++++++++++++++++++++++++++++++++\n" << std::endl;
|
|
||||||
std::cerr << "WARNING (ignored): Cannot do Schur complement with respect to non-diagonal block." << std::endl;
|
// Use sparse LU to solve the block submatrices i.e compute inv(D)
|
||||||
//OPM_THROW(std::logic_error, "Cannot do Schur complement with respect to non-diagonal block.");
|
const Eigen::SparseLU< M > solver(Jn[n]);
|
||||||
}
|
|
||||||
V diag = D.diagonal();
|
// compute inv(D)*bn for the update of the right hand side
|
||||||
Eigen::DiagonalMatrix<double, Eigen::Dynamic> invD = (1.0 / diag).matrix().asDiagonal();
|
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) {
|
for (int eq = 0; eq < num_eq; ++eq) {
|
||||||
if (eq == n) {
|
if (eq == n) {
|
||||||
continue;
|
continue;
|
||||||
}
|
}
|
||||||
|
const std::vector<M>& Je = eqs[eq].derivative();
|
||||||
|
const M& B = Je[n];
|
||||||
|
|
||||||
jacs[eq].reserve(num_eq - 1);
|
jacs[eq].reserve(num_eq - 1);
|
||||||
const M& B = eqs[eq].derivative()[n];
|
|
||||||
for (int var = 0; var < num_eq; ++var) {
|
for (int var = 0; var < num_eq; ++var) {
|
||||||
if (var == n) {
|
if (var == n) {
|
||||||
continue;
|
continue;
|
||||||
}
|
}
|
||||||
|
|
||||||
// Create new jacobians.
|
// Create new jacobians.
|
||||||
M schur_jac = eqs[eq].derivative()[var] - B * (invD * eqs[n].derivative()[var]);
|
// Add A
|
||||||
jacs[eq].push_back(schur_jac);
|
jacs[eq].push_back(Je[var]);
|
||||||
|
M& J = jacs[eq].back();
|
||||||
|
// solve Du = C
|
||||||
|
const M& u = solver.solve(Jn[var]);
|
||||||
|
// Subtract Bu (B*inv(D))
|
||||||
|
J -= B * u;
|
||||||
}
|
}
|
||||||
// Update right hand side.
|
// Update right hand side.
|
||||||
vals[eq] = eqs[eq].value().matrix() - B * (invD * eqs[n].value().matrix());
|
vals[eq] = eqs[eq].value().matrix() - B * Dibn;
|
||||||
}
|
}
|
||||||
|
|
||||||
// Create return value.
|
// Create return value.
|
||||||
@ -292,31 +302,26 @@ namespace Opm
|
|||||||
{
|
{
|
||||||
// The equation to solve for the unknown y (to be recovered) is
|
// The equation to solve for the unknown y (to be recovered) is
|
||||||
// Cx + Dy = b
|
// Cx + Dy = b
|
||||||
// y = inv(D) (b - Cx)
|
// Dy = (b - Cx)
|
||||||
// where D is the eliminated block, C is the jacobian of
|
// where D is the eliminated block, C is the jacobian of
|
||||||
// the eliminated equation with respect to the
|
// the eliminated equation with respect to the
|
||||||
// non-eliminated unknowms, b is the right-hand side of
|
// non-eliminated unknowms, b is the right-hand side of
|
||||||
// the eliminated equation, and x is the partial solution
|
// the eliminated equation, and x is the partial solution
|
||||||
// of the non-eliminated unknowns.
|
// of the non-eliminated unknowns.
|
||||||
// We require that D is diagonal.
|
|
||||||
|
|
||||||
// Find inv(D).
|
|
||||||
const M& D = equation.derivative()[n];
|
const M& D = equation.derivative()[n];
|
||||||
if (!isDiagonal(D)) {
|
|
||||||
std::cerr << "WARNING (ignored): Cannot do Schur complement with respect to non-diagonal block." << std::endl;
|
|
||||||
//OPM_THROW(std::logic_error, "Cannot do Schur complement with respect to non-diagonal block.");
|
|
||||||
}
|
|
||||||
V diag = D.diagonal();
|
|
||||||
Eigen::DiagonalMatrix<double, Eigen::Dynamic> invD = (1.0 / diag).matrix().asDiagonal();
|
|
||||||
|
|
||||||
// Build C.
|
// Build C.
|
||||||
std::vector<M> C_jacs = equation.derivative();
|
std::vector<M> C_jacs = equation.derivative();
|
||||||
C_jacs.erase(C_jacs.begin() + n);
|
C_jacs.erase(C_jacs.begin() + n);
|
||||||
ADB eq_coll = collapseJacs(ADB::function(equation.value(), C_jacs));
|
ADB eq_coll = collapseJacs(ADB::function(equation.value(), C_jacs));
|
||||||
const M& C = eq_coll.derivative()[0];
|
const M& C = eq_coll.derivative()[0];
|
||||||
|
|
||||||
|
// Use sparse LU to solve the block submatrices
|
||||||
|
const Eigen::SparseLU< M > solver(D);
|
||||||
|
|
||||||
// Compute value of eliminated variable.
|
// Compute value of eliminated variable.
|
||||||
V elim_var = invD * (equation.value().matrix() - C * partial_solution.matrix());
|
const Eigen::VectorXd b = (equation.value().matrix() - C * partial_solution.matrix());
|
||||||
|
const Eigen::VectorXd elim_var = solver.solve(b);
|
||||||
|
|
||||||
// Find the relevant sizes to use when reconstructing the full solution.
|
// Find the relevant sizes to use when reconstructing the full solution.
|
||||||
const int nelim = equation.size();
|
const int nelim = equation.size();
|
||||||
|
Loading…
Reference in New Issue
Block a user