is pc = 0? Otherwise I don’t understand the boundary condition. Your permeability seems to be a matrix. I that correctly implemented? Should be a diagonal matrix with the identical entries according to your implementation.
If the numeric differentiation does not work, I would rather suspect that something with your setup is wrong. Non-physical boundary / initial conditions or similar. 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 31. Mar 2020, at 19:35, Dmitry Pavlov <[email protected]> wrote: > > >> >>> >>> >>> 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 [email protected] https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
