Physical MPCs#

Coupling is an art, and often requires special off-diagonal work for globally implicit coupling, and fancy games can be played with domains to couple across domain interfaces both implicitly and sequentially. Physical MPCs derive from default MPCs to provide special implementations of some methods.

Coupled Water MPC#

A coupler which integrates surface and subsurface flow implicitly.

Couples Richards equation to surface water through continuity of both pressure and fluxes. This leverages subsurface discretizations that include face-based unknowns, and notes that those face unknowns that correspond to surface faces are co-located with the surface cell pressure, and therefore are equivalent. In this approach (described in detail in a paper that is in review), the surface equations are directly assembled into the subsurface discrete operator.

mpc-coupled-water-spec

  • “PKs order” [Array(string)] The use supplies the names of the coupled PKs. The order must be {subsurface_flow_pk, surface_flow_pk} (subsurface first).

  • “subsurface domain name” [string] domain

  • “surface domain name” [string] surface

  • “water delegate” [mpc-delegate-water-spec] A `Coupled Water Globalization Delegate`_ spec.

INCLUDES:

Coupled Cells MPC#

A coupler which solves two PDEs on the same domain.

This is a StrongMPC which uses a preconditioner in which the block-diagonal cell-local matrix is dense. If the system looks something like:

A( y1, y2, x, t ) = 0 B( y1, y2, x, t ) = 0

where y1,y2 are spatially varying unknowns that are discretized using the MFD method (and therefore have both cell and face unknowns), an approximation to the Jacobian is written as

[ dA_c/dy1_c dA_c/dy1_f dA_c/dy2_c 0 ] [ dA_f/dy1_c dA_f/dy1_f 0 0 ] [ dB_c/dy1_c 0 dB_c/dy2_c dB_c/dy2_f ] [ 0 0 dB_f/dy2_c dB_f/dy2_f ]

Note that the upper left block is the standard preconditioner for the A system, and the lower right block is the standard precon for the B system, and we have simply added cell-based couplings, dA_c/dy2_c and dB_c/dy1_c.

Most commonly this is used to couple flow and energy equations on the same mesh. In the temperature/pressure system, these extra blocks correspond to

\[\frac{\partial \Theta}{\partial T} \; , \; \frac{\partial E}{\partial p}\]

mpc-coupled-cells-spec

  • “domain name” [string] Domain of simulation

  • “conserved quantity A” [string] Key of the first sub-PK’s conserved quantity.

  • “conserved quantity B” [string] Key of the second sub-PK’s conserved quantity.

  • “primary variable A” [string] Key of the first sub-PK’s primary variable.

  • “primary variable B” [string] Key of the second sub-PK’s primary variable.

  • “no dA/dy2 block” [bool] false Excludes the dA_c/dy2_c block above.

  • “no dB/dy1 block” [bool] false Excludes the dB_c/dy1_c block above.

INCLUDES:

Subsurface MPC#

A coupler which solves flow and energy in the subsurface.

This MPC provides most nearly all terms for an approximate Jacobian for coupling three-phase Richards equation (the `Permafrost Flow PK`_) to the three-phase Energy equation (the `Three-Phase subsurface Energy PK`_).

Many options are provided for turning on and off various aspects of this Jacobian, so it is useful to mathematically write out these terms. The equations are:

\[\begin{split}\frac{\partial \Theta}{\partial t} - \nabla \frac{k_r n_l}{\mu} K ( \nabla p + \rho g \cdot \hat{z} ) = Q_w \\ \frac{\partial E}{\partial t} - \nabla \cdot \kappa \nabla T + \nabla \cdot \mathbf{q} e(T) = Q_w e(T) + Q_e\end{split}\]

Note that all of the following are dependent on \(p\) and/or \(T\):

\[\Theta(p,T), k_r(p,T), n_l(p,T), \mu(T), \rho(p,T), E(p,T), \kappa(p,T), e(T)\]

Also, both source terms \(Q_w\) and \(Q_e\) may or may not depend on \(p\) and \(T\).

Note also that the Darcy flux \(\mathbf{q}\) used in the advection of energy is given by the Darcy flux:

\[\mathbf{q} = -\frac{k_r n_l}{\mu} K ( \nabla p + \rho g \cdot \rho g \cdot \hat{z} )\]

Differentiating these two equations in their two unknowns gives the following four blocks in the approximate Jacobian:

\(\frac{\partial F_1}{\partial p}\): this is the Richards equation diagonal block, and is controlled inside that PK.

\(\frac{\partial F_1}{\partial T}\) includes terms for:

  • \(\frac{\partial \Theta}{\partial T}\) This term is the cell-local diagonal block.

  • The partial derivative of the divergence of the Darcy flux with respect to temperature is dominated by \(\frac{\partial}{\partial T} \frac{k_r n_l}{\mu}\). This is because the relative permeability is strongly dependent upon phase change (the freezing equals drying approximation). This term is referred to as the “d div q / dT” term.

\(\frac{\partial F_2}{\partial p}\) includes terms for:

  • \(\frac{\partial E}{\partial p}\) This term is the cell-local diagonal block.

  • The partial derivative of the energy diffusion term with respect to pressure is dominated by \(\frac{\partial \kappa}{\partial p}\) through phase change – at a constant temperature, but changing pressure, phase change can result in large changes to thermal conductivity. This is referred to as the “div K grad T / dp” term.

\(\frac{\partial F_2}{\partial T}\): this is the energy equation diagonal block, and is controlled inside that PK.

Also, at this level, where we know more about the flux used in the energy equation (it is the Darcy flux), we can do a better approximation of the derivative of the advection of energy term with respect to both temperature and pressure. For instance, enthalpy is only weakly dependent on pressure, so we can use the derivative of the divergence of the Darcy flux with respect to pressure (from the Richards block) in the advection term in the \(\frac{\partial F_2}{\partial p}\) block, and approximate \(\frac{\partial k_r}{\partial T}\) in the advection term as well. These terms are referred to as “div hq / dp,T terms”. Note the missing initial “d” here relative to other terms.

The behavior of this MPC’s preconditioner can be set by an option, “preconditioner type”. Really users should not change this from the default, except in expert cases or for comparison’s sake, but the options are:

  • “picard” is the default, this uses all available terms, and enables the “suppress” options for finer-grained control.

  • “none” No preconditioner never works.

  • “block diagonal” This is what one would get from the default StrongMPC_. This probably never works.

  • “no flow coupling” This keeps the accumulation terms, but turns off all the non-local blocks. This is equivalent to Coupled Cells MPC.

  • “ewc” CURRENTLY DEPRECATED/BROKEN/DISABLED In addition to the “picard” coupling, this also always does a change of variables, whereby we first invert to calculate primary variable corrections, then do a change of variables to calculate the linearized corrections in energy and water content space. We then apply those corrections, and invert to find the primary variable changes that would have made those corrections. This is called the “energy and water content” algorithm, and is related to similar variable changing approaches by Krabbenhoft (for flow) and Knoll (for energy), but in the multivariate approach. This is somewhat bad, becuase while it fixes some corrections, it breaks others.

  • “smart ewc” CURRENTLY DEPRECATED/BROKEN/DISABLED Does the “ewc” algorithm above, but tries to be smart about when to do it. This algorithm helps when we are about to fall off of the latent heat cliff. If we can guess when to do it, we have a better chance of not breaking things. This seems like it ought to be helpful, but often doesn’t do as much as one might hope.

Note this “ewc” algorithm is just as valid, and more useful, in the predictor (where it is not deprecated/disabled). There, we extrapolate a change in pressure and temperature, but often do better to extrapolate in water content and energy space, then invert (locally) for pressure and temperature corrections that meet that extrapolation. Both of these globalization algorithms are supported by the `EWC Globalization Delegate`_ object.

mpc-subsurface-spec

  • “domain name” [string] Domain of simulation

  • “preconditioner type” [string] picard See the above for detailed descriptions of the choices. One of: “none”, “block diagonal”, “no flow coupling”, “picard”, “ewc”, and “smart ewc”.

  • “supress Jacobian terms: div hq / dp,T” [bool] false If using picard or ewc, do not include this block in the preconditioner.

  • “supress Jacobian terms: d div q / dT” [bool] false If using picard or ewc, do not include this block in the preconditioner.

  • “supress Jacobian terms: d div K grad T / dp” [bool] false If using picard or ewc, do not include this block in the preconditioner.

  • “ewc delegate” [mpc-delegate-ewc-spec] A `EWC Globalization Delegate`_ spec.

INCLUDES:

Surface MPC#

Permafrost MPC#

A coupler which solves flow and energy both surface and subsurface.

This MPC handles the coupling of surface energy and flow to subsurface energy and flow for integrated hydrology with freeze/thaw processes.

mpc-permafrost-spec

  • “PKs order” [Array(string)] The user supplies the names of the coupled PKs. The order must be {subsurface_flow_pk, subsurface_energy_pk, surface_flow_pk, surface_energy_pk}.

  • “subsurface domain name” [string] domain

  • “surface domain name” [string] surface

  • “mass exchange flux key” [string] SURFACE_DOMAIN-surface_subsurface_flux

  • “energy exchange flux key” [string] SURFACE_DOMAIN-surface_subsurface_energy_flux

  • “water delegate” [mpc-delegate-water-spec] A `Coupled Water Globalization Delegate`_ spec.

INCLUDES: