Dear Denis,
I suggest to continue the issue #168 in this topic,
since the original issue seems to be more fundamental and is not really related to just Cutheel-mckee algorithm.
I am implementing a scheme of numerical solution of the two-phase oil-reservoir model. The phases are oil and water.
As a test case I consider the simplest system, both fluids are incompressible, densities and viscosities are fixed. No capillary pressure.
I am using the Finite Volume Method for the partial differential equations discretization, and
I am solving the problem in (s_water; pressure) variables using the Newton algorithm to linearize the algebraic problem.
s_oil is simply s_oil = 1 - s_water.
The particular problem that I report here is a 38x38 cells 2D square reservoir with two wells.
The initial condition is s_oil_0 = 1, p_0 = -1. The boundary conditions are left and right sides are impermeable for liquids,
the upper and lower bounds of the box are at fixed pressure, p = p_0 = -1. One well is is injection well and the other is a production well.
The vector of unknows is composed as s_water for every cell followed by pressure for every cell.
Thus a matrix with a 2x2 block structure is obtained as shown in figure
The lower-right block contains some non-zero elements due to boundary condition which currently assume an inflow of water through the upper and lower boundaries of reservoir.
The problem is that while using the AMGCL library, the linear system can be solved for the grids with the number of cells less than 39 in every direction. The finest allowed grid is of 38x38 cells, and the lower right block of the matrix contains 36 successive lines of zeros between the lines corresponding to the boundary conditions (with non-zero elements on the main diagonal).
I am using the Eigen library and the RowMajor SparseMatrix class to store my matrix A
typedef Eigen::SparseMatrix<double, RowMajor> SpMat;
SpMat A(2*N, 2*N);
Then, I use the following code to interface with the AMGCL
#include <amgcl/adapter/eigen.hpp>
#include <amgcl/backend/builtin.hpp>
#include <amgcl/make_solver.hpp>
#include <amgcl/solver/gmres.hpp>
#include <amgcl/amg.hpp>
#include <amgcl/coarsening/smoothed_aggregation.hpp>
#include <amgcl/relaxation/spai0.hpp>
typedef amgcl::make_solver<
amgcl::amg<
amgcl::backend::builtin<double>,
amgcl::coarsening::smoothed_aggregation,
amgcl::relaxation::spai0
>,
amgcl::solver::gmres<amgcl::backend::builtin<double> >
> Solver_AMG;
Solver_AMG solve(A);
int iters;
double error;
std::tie(iters, error) = solve(b0, u_Next);
For the number of grid cells more or equal to 39x39 I obtain the following error
Unhandled exception at 0x00007FFE31883E49 in OilReservoir_Eigen.exe: Microsoft C++ exception: std::runtime_error at memory location 0x0000002F1D0FE908. occurred
Pressing F5 in Visual Studio 2019, I get this with D of all zeros (715 elements total)
OilReservoir_Eigen.exe has triggered a breakpoint. occurred
So, my concerns are as follows.
I appreciate if you could kindly comment on them.
-
The Eigen library implements SparseLU direct solver that solves the problem for any grid sizes.
And the solution is physically reasonable.
So, I do not prone to believe that the issue is "This looks like a memory error somewhere in your program",
as you suggested in issue #168.
-
I do not think that the Schur pressure correction approach is (physically) justified approach in this case.
Consider the reservoir with s_oil = s_water = 0.5 in every cell and uniform pressure. Then, the two lines of blocks would be
equivalent, so I cannot simply disregard the non-zeros in the lower right block of matrix in the iterative solver.
-
Should I reorder the vector of unknows so that every two successive lines in the matrix are related to
the same grid cell, but represent two equations, one for every phase. Then I could say that the matrix contains NxN cells (N - totl number of cells in the grid), where every cell is a 2x2 block. Thus, every grid cell is described by a 2x2 "tensor". Then, I could use the option amgcl::static_matrix<2,2,double> to initialize the amg solver?
-
Does the implemented backend for Eigen support the Eigen::SparseMatrix<double, RowMajor> simultaneously with amgcl::static_matrix<2,2,double>?
2.5k Aug 5, 2022
33 Jul 14, 2022
59 Jul 18, 2022
692 Jul 28, 2022
35 Jul 5, 2022
12 Jul 28, 2022
571 Jul 31, 2022
66 Jul 24, 2022
16 Apr 19, 2022
491 Aug 5, 2022
184 Jun 17, 2022
600 Jul 29, 2022
728 Jul 26, 2022
523 Jul 30, 2022
4.7k Aug 1, 2022
3.4k Aug 9, 2022
385 Jul 26, 2022
879 Aug 5, 2022
18 May 2, 2022