Hi Dmitry,

have a look at 
https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/merge_requests/1952 
<https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/merge_requests/1952>
It should fix the bug for tpfa and also another bug in the box method 
implementation of the analytic Jacobian.
It also adds tests for the analytic Jacobian that pass against reference 
solutions created with numeric Jacobian.

Thanks a lot for helping us improve DuMux! This is highly appreciated.

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 6. Apr 2020, at 11:35, Dmitry Pavlov <[email protected]> wrote:
> 
> Dear Timo,
> 
> I actually had 1x4 grid, so the "inner" two elements were connected to walls 
> but not to inlet/outlet boundaries.
> 
> Here is a similar example (attached), but with 4x4 grid, with capillary 
> pressure restored, and with initial sn = 0.5. Hope it helps.
> Best regards,
> 
> Dmitry
> 
> 
> On 05.04.2020 12:40, Timo Koch wrote:
>> 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] 
>> <mailto:[email protected]>
>> D-70569 Stuttgart             url: www.iws.uni-stuttgart.de/en/lh2/ 
>> <http://www.iws.uni-stuttgart.de/en/lh2/>
>> _________________________________________________
>> 
>>> On 4. Apr 2020, at 12:48, Dmitry Pavlov <[email protected] 
>>> <mailto:[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] <mailto:[email protected]>
>>> https://listserv.uni-stuttgart.de/mailman/listinfo/dumux 
>>> <https://listserv.uni-stuttgart.de/mailman/listinfo/dumux>
>> 
> <4x4-analytic.txt><4x4-numeric.txt>_______________________________________________
> 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

Reply via email to