Timo,
My analytic Jacobian may or may not be correct, I agree. But it is
fairly close to the numeric one and the deltaU is very similar.
Now, I really can not understand what is happening with this phase
composition system. I assumed that whatever problems I have in the
Jacobian, they will manifest themselves when I solve the linear system
with the Jacobian and not when I solve some other system. I do not see
anything non-physical in pressure and saturation values (see below).
More to say, I do not have temperature equations at all. Component
densities do not change. Fluids are not compressible. I do not have
phase transition. There is no gas nor solid phase. Even the diffusion is
turned off. Though I am using 2pnc model, I do not use all functionality
that it provides. Not in this example, anyway. Finally, the third
component (soluble) is not actually present in this example. So I can
not imagine what can make the system non-physical.
Here is debug output that I printed out from
MiscibleMultiPhaseComposition::solve before error happened:
pressure(0) = 9.313226e-10
saturation(0) = 2.121000e-01
pressure(1) = 1.200000e+06
saturation(1) = 7.879000e-01
M =
0 0 0 -1.2e+08 0 0
0 9.31323e-08 0 0 -0 0
0 0 0 0 0 -1.2e+08
1 1 1 0 0 0
0 0 0 1 1 1
0 0 1 0 0 0
M after preconditioning =
0 0 0 -12000 0 0
0 9.31323e-12 0 0 -0 0
0 0 0 0 0 -1.2
1 1 1 0 0 0
0 0 0 1 1 1
0 0 1 0 0 0
b = 0 0 0 1 1 0
x after failed solution = 1 0 0 0 1 0
The matrix is not singular, though it is ill-conditioned. It is
interesting that FieldMatrix::solve() does find the correct solution
before throwing the Dune::FMatrixError.
Now, about those numbers -1.2e+08 and 9.31323e-08 whose ratio probably
causes the solve() to fail. I discovered that they come from the
fluidState.fugacityCoefficient() call. And I am not sure if I the
concept of fugacity is important for my program at all, given the
underlying physical model.
Here is what I once wrote in my custom implementation of the fluid system.
template <class FluidState>
static Scalar fugacityCoefficient(const FluidState &fluidState,
int phaseIdx,
int compIdx)
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
assert(0 <= compIdx && compIdx < numComponents);
if (phaseIdx == phase0Idx)
if (compIdx == wettingComp0Idx || compIdx == wettingComp1Idx)
return 0.0; // water and soluble stay in wetting phase
if (phaseIdx == phase1Idx)
if (compIdx == nonwettingCompIdx)
return 0.0; // oil stays in oil phase
return 100.0; // arbitrary?
}
You see that I am not sure about what I should return if phaseIdx is
nonwetting phase and compIdx is water/soluble, or if phaseIdx is wetting
phase and compIdx is oil. Maybe this causes all the trouble with
MiscibleMultiPhaseComposition::solve() ?
On 27.05.2020 15:32, Timo Koch wrote:
Hi Dmitry,
the error is from the constraint solver that solves for the phase
composition given some state variables (e.g. pressure, temperature).
This is called in the volume variables.
Most likely there is something wrong in your analytic Jacobian which
leads to a non-physical solution and then the phase composition cannot
be determined.
The constraint solver solves a small local linear equation system,
this is I guess why the error message say “Matrix”. But the error
message is a bit unfortunate since we don’t solve matrices but linear
equation systems.
I’ll have a look and maybe propose a better message.
Timo
On 26. May 2020, at 22:47, Dmitry Pavlov <[email protected]
<mailto:[email protected]>> wrote:
Hello,
I am a bit puzzled about this error:
[solve:/home/dpavlov/DUMUX/dumux/dumux/material/constraintsolvers/misciblemultiphasecomposition.hh:182]:
Matrix for composition of phases could not be solved.
I have two phase (wetting + nonwetting), three component (water + oil
+ soluble) problem. I am using the box method. The error arises even
if the concentration of the soluble is zero in both initial and
boundary conditions. It arises with non-zero concentration as well. I
am not sure I understand what it means. It happens with analytic
Jacobian but not with numeric, though the two Jacobians are close and
the resulting deltaU are close.
[...]
_______________________________________________
Dumux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
