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
