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.
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?
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?
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.
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;
}
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
