Dear Timo,
addRobinFluxDerivatives() is only needed if you want to implement an
exact / analytical Jacobian matrix for your problem. In the default
setting DuMux computes the Jacobian using numerical differentiation
and doesn’t need any user-implemented derivatives.
I am having trouble now with Newton method not converging sometimes. The
convergence depends on the grid size and on the influx input parameter,
and I am suspecting the numerical differentiation may be the cause. I
played with Assembly.NumericDifference.BaseEpsilon parameter, to no avail.
So I decided to check out the analytical addRobinFluxDerivatives(). I
created a MyLocalResidual class, inherited from
TwoPIncompressibleLocalResidual, and implemented
addRobinFluxDerivatives() method in it. I also use Assembler =
FVAssembler<TypeTag, DiffMethod::analytic> instead of
DiffMethod::numeric, of course.
The program compiles, and does call my method, but the solution diverges
dramatically. I checked the derivatives against ones that I calculated
numerically, and they are fine, to my understanding anyway. I suspect
that maybe I miss something about the units of the derivatives that
addRobinFluxDerivatives is supposed to add.
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 code snippet is given below for convenience.
template<class PartialDerivativeMatrices>
void addRobinFluxDerivatives(PartialDerivativeMatrices&
derivativeMatrices,
const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables&
curElemVolVars,
const ElementFluxVariablesCache&
elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
const auto insideScvIdx = scvf.insideScvIdx();
const auto& globalPos = scvf.ipGlobal();
auto ccGlobalPos = element.geometry().center();
auto& dI_dI = derivativeMatrices[insideScvIdx];
if (globalPos[0] > problem.bBoxMax()[0] - Problem::eps_) {
const MaterialLawParams& mlp =
problem.spatialParams().materialLawParamsAtPos(ccGlobalPos);
const auto& volVars = curElemVolVars[scvf.insideScvIdx()];
Scalar pw = volVars.pressure(0);
Scalar sw = volVars.saturation(0);
Scalar krw = MaterialLaw::krw(mlp, sw);
Scalar krn = MaterialLaw::krn(mlp, sw);
const auto& insideScv = fvGeometry.scv(insideScvIdx);
auto dKrw_dSn = -1.0*MaterialLaw::dkrw_dsw(mlp, sw);
auto dKrn_dSn = -1.0*MaterialLaw::dkrn_dsw(mlp, sw);
Scalar gradp = (problem.pwAtBoundary_ - pw) / (globalPos[0] -
ccGlobalPos[0]);
Scalar dgradp_dpw = - 1.0 / (globalPos[0] - ccGlobalPos[0]);
FieldMatrix K =
problem.spatialParams().permeabilityAtPos(ccGlobalPos);
dI_dI[0][Indices::pressureIdx] += - volVars.density(0) *
volVars.mobility(0) * K[0][0] * dgradp_dpw;
dI_dI[0][Indices::saturationIdx] += - volVars.density(0) *
dKrw_dSn / volVars.fluidState().viscosity(0) * K[0][0] * gradp;
dI_dI[1][Indices::pressureIdx] += - volVars.density(1) *
volVars.mobility(1) * K[0][0] * dgradp_dpw;
dI_dI[1][Indices::saturationIdx] += - volVars.density(1) *
dKrn_dSn / volVars.fluidState().viscosity(1) * K[0][0] * gradp;
}
}
The corresponding code snippet from neumann () is:
const auto& volVars = elemVolVars[scvf.insideScvIdx()];
auto ccGlobalPos = element.geometry().center();
FieldMatrix K =
this->spatialParams().permeabilityAtPos(ccGlobalPos);
Scalar pw = volVars.pressure(0);
Scalar dp_dn = (pwAtBoundary_ - pw) / (globalPos[0] -
ccGlobalPos[0]);
values[0] = - volVars.density(0) * volVars.mobility(0) *
K[0][0] * dp_dn;
values[1] = - volVars.density(1) * volVars.mobility(1) *
K[0][0] * dp_dn;
Best regards,
Dmitry
_______________________________________________
Dumux mailing list
Dumux@listserv.uni-stuttgart.de
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
