Merge branch 'fix/dkrn-dswe-brookscorey' into 'master'

Fix wrong implementation of dkrn_dswe in Brooks-Corey and Van Genuchten

Closes #844

See merge request !1912

CHANGELOG.md

@@ -20,6 +20,7 @@ An additional new option is `Vtk.CoordPrecision` which changes the precision of
 - Remove `Grid.HeapSize` as dune-ugrid removed the according feature as well.
 - __Van Genuchten__: Corrected VanGenuchten-Mualem exponent in the non-wetting saturation formula (`1/3` instead of `1/2` (or `l`, see above))
 - __Van Genuchten__: Corrected VanGenuchten-Mualem implementation of `dkrn/dSw`
+- __Brooks-Corey__: Corrected Brooks-Corey implementation of `dkrn/dSw` and added the derivatives for the regularized version
 - __AMGBackend__: The internal structure of the AMGBackend and the ParallelISTLHelper has been overhauled, as only used by the AMG, we did not make the changes backwards-compatible
 - The global default parameters for linear solvers have been removed and moved to the class `LinearSolver`.
 This only affects users that directly obtain this parameter via `getParam` somewhere in the code.

dumux/material/fluidmatrixinteractions/2p/brookscorey.hh

@@ -269,10 +269,10 @@ public:
 
         swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
 
-        return 2.0*(swe - 1)*(1 + pow(swe, 2.0/params.lambda())*(1.0/params.lambda() + 1.0/2
-                                                                 - swe*(1.0/params.lambda() + 1.0/2)
-                                                                )
-                             );
+        const auto lambdaInv = 1.0/params.lambda();
+        const auto swePow = pow(swe, 2*lambdaInv);
+        return 2*(swe - 1.0)*(1.0 + (0.5 + lambdaInv)*swePow - (1.5 + lambdaInv)*swePow*swe);
+
     }
 
 };

dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh

@@ -280,7 +280,7 @@ public:
 
         const auto sne = 1.0 - swe;
         const auto x = 1.0 - pow(swe, 1.0/params.vgm());
-        return -pow(sne, params.vgl()-1.0) * pow(x, 2*params.vgm() - 1.0) * ( params.vgl()*x - 2.0*sne/swe*(1.0 - x) );
+        return -pow(sne, params.vgl()-1.0) * pow(x, 2*params.vgm() - 1.0) * ( params.vgl()*x + 2.0*sne/swe*(1.0 - x) );
     }
 
 };