Dear Dmitry, Thank you so much for your detailed investigation.
I’m pretty convinced from your output that the analytical Jacobian is currently wrong. (It is also currently untested.) In case you don’t know how to interpret the output: The matrix is blocked into primary variable blocks (in this case 2x2 sub matrices). So e.g. the first and second entries in row 0 and 1 are the subblock of the derivatives of element 0 residuals (wetting, non-wetting phase balance) with respect to the dofs in element 0. (0,0) -> drw/dpw (1,0) -> drn/dpw (0,1) -> drw/dsn (1,1) -> drn/dsn Looking at the storage derivative (twopincompressiblelocalresidual.hh) I’m pretty sure the sign is wrong for drn/dsn. There may be further mistakes in the fluxes. Your grid was 2x2 I assume from the matrix. That means that all element are connected to a boundary. Could you post the same output for a 4x4 grid? So we can better isolate how entries of inner cells differ? Also some effects don’t show if one of the saturations is zero because some term drop out. Can you start with sn=0.5? Best wishes Timo -- _________________________________________________ Timo Koch phone: +49 711 685 64676 IWS, Universität Stuttgart fax: +49 711 685 60430 Pfaffenwaldring 61 email: [email protected] D-70569 Stuttgart url: www.iws.uni-stuttgart.de/en/lh2/ _________________________________________________ > On 4. Apr 2020, at 12:48, Dmitry Pavlov <[email protected]> wrote: > > Hello, > > I just ran into a weird problem with DuMux 3.1. As I reported before, I was > having trouble replacing the numeric derivatives with my own analytic ones. > Now I made a very small and simplified example and made some debug printing > to track down the issue. > > The problem being solved is 2p (TwoPIncompressibleLocalResidual), fully > implicit scheme. The assembler is FVAssembler<TypeTag, DiffMethod::numeric>. > The boundary conditions are: neumann (injection well), all Dirichlet > (production well). No customary flux derivatives, no other fancy stuff. > > Here is the first output of the debug printing that I added into > NewtonSolver::solveLinearSystem_: > > Jacobian [n=4,m=4,rowdim=8,coldim=8] > row 0 5.20412e-12 -7.65000e+09 -5.10436e-12 0.00000e+00 . > . . . > row 1 2.43297e-09 6.50250e+09 -2.43297e-09 0.00000e+00 . > . . . > row 2 -5.10436e-12 0.00000e+00 1.02087e-11 -7.65000e+09 -5.10436e-12 > 0.00000e+00 . . > row 3 -2.43297e-09 0.00000e+00 4.86593e-09 6.50250e+09 -2.43297e-09 > 0.00000e+00 . . > row 4 . . -5.10436e-12 0.00000e+00 1.02087e-11 > -7.65000e+09 -5.10436e-12 0.00000e+00 > row 5 . . -2.43297e-09 0.00000e+00 4.86593e-09 > 6.50250e+09 -2.43297e-09 0.00000e+00 > row 6 . . . . -5.10436e-12 > 0.00000e+00 1.53131e-11 -7.65000e+09 > row 7 . . . . -2.43297e-09 > 0.00000e+00 7.29888e-09 6.50250e+09 > residual [blocks=4,dimension=8] > row 0 -2.17e-02 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 1.02e-06 > 4.87e-04 > > And the simulation goes fine after that (more or less). And here is what the > first debug output becomes when I use FVAssembler<TypeTag, > DiffMethod::analytic> instead: > > Jacobian [n=4,m=4,rowdim=8,coldim=8] > row 0 5.10435e-12 -7.65000e+09 -5.10435e-12 0.00000e+00 . > . . . > row 1 2.43296e-09 -6.50250e+09 -2.43296e-09 0.00000e+00 . > . . . > row 2 -5.10435e-12 0.00000e+00 1.02087e-11 -7.65000e+09 -5.10435e-12 > 0.00000e+00 . . > row 3 -2.43296e-09 0.00000e+00 4.86593e-09 -6.50250e+09 -2.43296e-09 > 0.00000e+00 . . > row 4 . . -5.10435e-12 0.00000e+00 1.02087e-11 > -7.65000e+09 -5.10435e-12 0.00000e+00 > row 5 . . -2.43296e-09 0.00000e+00 4.86593e-09 > -6.50250e+09 -2.43296e-09 0.00000e+00 > row 6 . . . . -5.10435e-12 > 0.00000e+00 1.53131e-11 -7.65000e+09 > row 7 . . . . -2.43296e-09 > 0.00000e+00 7.29889e-09 -6.50250e+09 > residual [blocks=4,dimension=8] > row 0 -2.17e-02 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 1.02e-06 > 4.87e-04 > > And the simulation stucks instantly. You see that the odd diagonal elements > have the opposite signs as compared to the ones of the analytic Jacobian. > > Did I run into a bug? Did I misuse DuMux? I will appreciate any comments. > > Best regards, > > Dmitry > > > _______________________________________________ > Dumux mailing list > [email protected] > https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
_______________________________________________ Dumux mailing list [email protected] https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
