I supposed that I should calculate the derivatives of the fluxes (in
kg / m^2 * s, as they are calculated by neumann()) w. r. t. the
wetting phase pressure (Pa) and nonwetting phase saturation. Am I
correct?
The fluxes are integrated over the sub control volume face for finite
volumes so the unit should be kg/s, I think. You need to multiply with
the flux derivatives with the face area (scvf.area()).
Oh, OK. Thanks you. I saw in the sources that neumann()'s output gets
multiplied by the face area inside DuMux itself, but was not sure about
the derivatives.
I just multiplied all the four derivatives by scvf.area(). It got a lot
better, but it still diverges, though not as dramatically as before. So
I kind of ran out of ideas now. Maybe I should multiply
addRobinFluxDerivatives's output by porosity? Or neumann's?
Here is the output accompanied by some debug printing. On the same
setting, an assembler DiffMethod::numeric does not diverge.
Assemble: r(x^k) = dS/dt + div F - q; M = grad r sw = 0.200000, pw =
100000.000000
Newton iteration 1 done, maximum relative shift = 1.762e+00
Assemble: r(x^k) = dS/dt + div F - q; M = grad r sw = 0.200000, pw =
106945.166613
Newton iteration 2 done, maximum relative shift = 5.337e+00
Assemble: r(x^k) = dS/dt + div F - q; M = grad r sw = 0.200000, pw =
-235123.442414
Newton iteration 3 done, maximum relative shift = 2.465e+03
Assemble: r(x^k) = dS/dt + div F - q; M = grad r sw = 0.200000, pw =
421258.254293
Newton iteration 4 done, maximum relative shift = 3.306e-01
Assemble: r(x^k) = dS/dt + div F - q; M = grad r sw = 0.200000, pw =
421245.329986
Newton iteration 5 done, maximum relative shift = 3.754e+00
Assemble: r(x^k) = dS/dt + div F - q; M = grad r sw = 0.200000, pw =
421234.078144
Newton iteration 6 done, maximum relative shift = 5.252e+00
Assemble: r(x^k) = dS/dt + div F - q; M = grad r sw = 0.200000, pw =
421211.216791
Newton iteration 7 done, maximum relative shift = 2.364e+00
Assemble: r(x^k) = dS/dt + div F - q; M = grad r sw = 0.200000, pw =
421239.307957
Newton iteration 8 done, maximum relative shift = 3.006e+00
Assemble: r(x^k) = dS/dt + div F - q; M = grad r sw = 0.200000, pw =
421233.214948
Newton iteration 9 done, maximum relative shift = 1.382e+00
Assemble: r(x^k) = dS/dt + div F - q; M = grad r sw = 0.200000, pw =
421232.752614
Newton iteration 10 done, maximum relative shift = 1.133e+01
Assemble: r(x^k) = dS/dt + div F - q; M = grad r sw = 0.200000, pw =
421197.730506
Newton iteration 11 done, maximum relative shift = 5.195e+00
Assemble: r(x^k) = dS/dt + div F - q; M = grad r sw = 0.200000, pw =
421149.236607
Newton iteration 12 done, maximum relative shift = 4.703e+00
Assemble: r(x^k) = dS/dt + div F - q; M = grad r sw = 0.200000, pw =
421205.601966
Newton iteration 13 done, maximum relative shift = 1.536e+02
Assemble: r(x^k) = dS/dt + div F - q; M = grad r sw = 0.200000, pw =
421160.385747
Newton iteration 14 done, maximum relative shift = 4.639e-01
Assemble: r(x^k) = dS/dt + div F - q; M = grad r sw = 0.200000, pw =
421289.418459
Newton iteration 15 done, maximum relative shift = 4.101e+00
Assemble: r(x^k) = dS/dt + div F - q; M = grad r sw = 0.200000, pw =
421236.003036
Newton iteration 16 done, maximum relative shift = 5.802e+00
Assemble: r(x^k) = dS/dt + div F - q; M = grad r sw = 0.200000, pw =
421123.366818
Newton iteration 17 done, maximum relative shift = 1.621e+00
Assemble: r(x^k) = dS/dt + div F - q; M = grad r sw = 0.200000, pw =
422616.393919
Newton iteration 18 done, maximum relative shift = 1.197e+00
Newton solver did not converge with dt = 0.001 seconds. Retrying with
time step of 0.0005 seconds
Also: which DuMux version are you using? We just fixed some error in
the implementation of dkrn_dsw today.
I use DuMux 3.1. I have my own implementation of dkrn_dsw, and with the
sw at 0.2, the derivative of krn and krw are zero.
Best regards,
Dmitry
_______________________________________________
Dumux mailing list
Dumux@listserv.uni-stuttgart.de
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
