Dear Wencheng, Timo has already mentioned all the important details. So, I can only add a few minor points. If flow occurs in the soil and the macropore domain, the model is called a dual-permeability model. In a dual-porosity model, only flow in the macropores is considered. Perhaps this explains the observed differences between your DuMux dual-permeability model and the HYDRUS-1D dual-porosity model. I had implemented a two-phase (water/air) dual-permeability model in DuMux in 2010 (this model does no longer exist in DuMux). Therefore I extended the equations and defined the soil parameters for each domain (as Timo mentioned). The domains where coupled via the source term, which worked very well. I also conducted a comparison between a 1D dual-permeability model based on the Richards equation and the two-phase water/air counterpart (dual-podual-pohttps://doi.org/10.1007/s11242-012-0064-3 http://dx.doi.org/10.14279/depositonce-4972 , http://dx.doi.org/10.14279/depositonce-4972). Perhaps the results of these studies could be used for validation or comparison purposes. Best regards Leo
> Koch, Timo <[email protected]> hat am 13.01.2026 10:38 CET > geschrieben: > > > Dear Wencheng, > > welcome to the mailing list! I’m glad you found DuMux for this task and find > it useful. > I’ll try to answer your questions below inline. > > > > On 8 Jan 2026, at 21:03, Qi, Wencheng <[email protected]> wrote: > > > > Dear DuMuX Developer Team, > > > > Happy New Year! Thank you for your continued maintenance and development of > > DuMuX. My name is Wencheng Qi, and I am currently a postdoctoral researcher > > at the University of Illinois Urbana–Champaign working on soil–plant > > interactions. I am developing a dual-porosity model using the DuMuX > > Multidomain framework, and I have some questions regarding Jacobian > > assembly and cross-domain derivatives in Multidomain. Thank you very much > > for your time, and I would greatly appreciate any guidance you could > > provide. > > > > The dual-porosity model I am trying to implement is conceptually similar to > > the root–soil two-domain coupling approach in dumux/test: the two continuum > > domains represent the matrix and macropore regions, and each domain solves > > a Richards equation with its own Problem. Unlike an explicit geometric > > partitioning (e.g., “fracture” vs “matrix”), my two domains are > > conceptually overlapping: at each spatial location, both matrix and > > macropore states exist, and mass is exchanged between the two domains via > > an exchange term. > > > > If the the domain are fully overlapping, so essentially you have a true > dual/multi continuum model, and your are also using the same discretization > method for all equations, it is usually more efficient and effective to > extend a given model (e.g. Richards) by more equations rather than realizing > this in a multi-domain framework. The reason is that the Jacobian derivatives > can be computed a bit more efficiently in this case. Also the solver can > utilize the block structure. However, this is not to say that it doesn’t work > with multi-domain. With multidomain you get the Jacobian in global > model-specific block structure which might also be useful to you. > > > > > > After the model development, I tried to conduct a model validation between > > DuMuX and HYDRUS-1 (with dual-porosity). However, there is no agreement in > > flow variables (e.g., water exchange flux) between the two models. Thus, I > > have two major questions: > > > > 1. After adding the exchange term as a source term to the governing > > equation in each Problem, will this contribution automatically be included > > in evalLocalResidual? In particular, when DuMuX builds the Jacobian by > > numerical differentiation, will the newly added source term be > > differentiated and assembled into the Jacobian automatically? In other > > words, in typical DuMuX workflows, do I not need to manually implement the > > Jacobian specifically for this added source term > > > > Yes this is automatic (using numeric differences). In case you haven’t > overloaded the evalCouplingResidual function, per default, the full residual > will be evaluated which in includes all terms. > > > > 1. Because the exchange term depends on both the pressure in the current > > domain (p_this) and the pressure in the other domain (p_other), the > > Jacobian should, in principle, contain cross-domain block entries. Will > > these cross-domain derivatives also be included automatically? Or is it > > recommended to implement the exchange via a specific coupling > > residual/interface (rather than directly adding it in each Problem’s > > source(...)), or to explicitly provide/supplement Jacobian contributions > > somewhere in the CouplingManager/Assembler? > > > > > > Yes, the cross-domains are also computed automatically based on what you > provide in “couplingStencil”. The couplingStencil has to return a list of all > degrees of freedom of the other domain that have an influence on the local > residual of the given element in this domain. Often, if you forget some > degrees od freedom, you will obtain only a rough approximation to the > Jacobian and this increases the number of Newton iteration but you can often > still get convergence. However, this is not guaranteed. Note that the dofs > you return depend on the discretization scheme. For BOX/CVFE in your case it > would typically be all dofs of the element, whereas for CC-TPFA/MPFA it would > be the element dof. In case your coupling stencil only contains the dof of > the same control volume for TPFA, only the source/storage derivatives will be > evaluated correctly, not the fluxes, but I assume your model only couples via > source terms. > > > > > > > > Supplementary details > > > > My implementation workflow is as follows: I created two Richards > > subproblems (problem-matrix and problem-macropore) and implemented a > > CouplingManager to handle the mapping and exchange computation between the > > two domains. In the CouplingManager, I define a one-to-one DoF mapping > > (each grid cell corresponds to the cell at the same location in the other > > domain), and provide an exchangeSource(domain, scv) interface. In that > > function, I use the mapping to find the DoF in the other domain, read the > > pressure variables from curSol, compute the exchange term in the form > > q_mass = rho * alpha_p * (p_other − p_this) (where alpha_p is computed from > > my spatial parameters/hydraulic conductivity/geometric factors, etc.), and > > return the mass source density q_mass. Then, in each Problem I implement > > source(...) and add couplingManager().exchangeSource(...) to the continuity > > equation. > > > > My current understanding is: once the exchange term is added via > > source(...), it becomes one of the contributions to the local residual > > together with the storage term and the flux term. Therefore, when > > evalLocalResidual evaluates the total residual, it should evaluate these > > contributions together (storage + flux + source). If the Jacobian is then > > built by numerical differentiation, DuMuXperturbs the solution vector and > > re-evaluates the same total residual; in that process, the exchange source > > is re-evaluated as well and thus its derivatives should be captured > > automatically in the numerically differentiated Jacobian. Under this > > understanding, I would not need to add a separate Jacobian implementation > > specifically for the exchange source term. > > > > Yes this is correct. > > > > > > A simplified version of my CouplingManager exchangeSource is: > > > > // CouplingManager: compute exchange using mapped dofs and current solutions > > Scalar exchangeSource(domainI, scv) const > > { > > dofThis = scv.dofIndex(); > > dofOther = mapThisToOther[dofThis]; > > pThis = curSol(domainI)[dofThis ][pIdxThis ]; > > pOther = curSol(other )[dofOther][pIdxOther]; > > return q_mass = rho * alpha_p * (pOther - pThis); > > } > > > > > > > > And the source implementation in problem-macropore is: > > > > PrimaryVariables source(const Element& element, > > const FVElementGeometry& fvGeometry, > > const ElementVolumeVariables& elemVolVars, > > const SubControlVolume& scv) const > > { > > PrimaryVariables src(0.0); > > const Scalar q_mass = > > couplingManager().exchangeSource(Dune::index_constant<1>{}, scv); > > src[Indices::conti0EqIdx] += q_mass; > > return src; > > } > > > > > > The source in DuMux should be volume-specific. This might be the issue here. > If you have the full mass here, you need to divide by the scv.volume() to > have a volume-specific term. > (The term gets integrated and thus multiplied by the volume again later in > the local residual. But commonly/historically volume-specific source have > been a bit more convenient, not here.) > > Best wishes and let us know if you have any follow-up questions. > Timo > > P.S. You can also reach out to us through the matrix channel. > > —— > __________________________________________________________ > > Jun.-Prof. Dr.-Ing. Timo Koch > > Department for Hydromechanics and Modelling of Hydrosystems (LH2) > Institute for Modeling Hydraulic and Environmental Systems (IWS) > University of Stuttgart, Germany > > Institute webpage: https://www.iws.uni-stuttgart.de/en/lh2/ > Personal webpage: https://timokoch.github.io/ > Phone: +49 711 685 60461 > ___________________________________________________________ > > > > > > > If anything in my description or understanding is incorrect, I would > > greatly appreciate your corrections. I am also happy to provide additional > > implementation details to help you diagnose the issue. > > > > Thank you very much! > > > > Best regards, > > Wencheng Qi, PhD > > Postdoc Research Scientist > > Department of Crop Sciences > > University of Illinois Urbana–Champaign > > _______________________________________________ > > DuMux mailing list > > [email protected] > > https://listserv.uni-stuttgart.de/mailman/listinfo/dumux > > > _______________________________________________ > DuMux mailing list > [email protected] > https://listserv.uni-stuttgart.de/mailman/listinfo/dumux > Im Auftrag Dr.-Ing. Leopold Stadler -- Referat Daten und Methoden Abteilung Wasserbau im Binnenbereich Bundesanstalt für Wasserbau Federal Waterways Engineering and Research Institute Kußmaulstraße 17 | 76187 Karlsruhe https://www.baw.de
_______________________________________________ DuMux mailing list [email protected] https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
