--
_________________________________________________
Timo Koch phone: +49 711 685 64676
IWS, Universität Stuttgart fax: +49 711 685 60430
Pfaffenwaldring 61 email: timo.k...@iws.uni-stuttgart.de
D-70569 Stuttgart url: www.iws.uni-stuttgart.de/en/lh2/
_________________________________________________
> On 31. Mar 2020, at 18:46, Dmitry Pavlov <dmitry.pav...@outlook.com> wrote:
>
> 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.
Hi Dmitry,
I wouldn’t suspect numeric differentiation to be the problem. But it could be…
>
> 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 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()).
Also: which DuMux version are you using? We just fixed some error in the
implementation of dkrn_dsw today.
Best wishes
Timo
>
> 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
_______________________________________________
Dumux mailing list
Dumux@listserv.uni-stuttgart.de
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
