-- _________________________________________________ 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