Multiscale and Multiresolution Methods

General Description

MrHyDE is designed to use heterogeneous computational architectures to enable concurrent multiscale modeling and simulation. The general form of a coupled multiscale system is:

$$ \frac{\partial u_c}{\partial t} + F_c(u_c,u_f) = 0, $$ $$ \frac{\partial u_f}{\partial t} + F_f(u_c,u_f) = 0, $$

where $u_c$ denotes the coarse-scale variables and $u_f$ denotes the fine scale variables. Either $F_c$ or $F_f$ may also depend on space and/or time, but we omit this explicit dependence for simplicity. If we assume that the fine scale variables can be eliminated (preferably locally), then we can solve an equation only involving the coarse scale variables:

$$ \frac{\partial u_c}{\partial t} + F_c(u_c,h(u_c)) = 0, $$

where the notation $h(u_c)$ represents the local nonlinear elimination of the fine scale information, which depends on the current coarse scale state. As we will soon see, this is related to the Schur complement in the case of steady-state linear operators. At this point, there are two possible interpretations of the coupled system and the reduced system. If the objective is to compute an approximation to 𝑢𝑐 that is informed by the fine scale information, then one interpretation is that this can be obtained by solving the reduced system. However, if the objective is to approximate the fine scale solution that is informed by global information, then we can interpret the reduced system as providing the time evolution of the coarse scale global coupling between fine scale models. The appropriate interpretation is problem dependent.

For simplicity, we will ignore any spatial discretizations throughout this section and proceed as if we are solving coupled systems of ODEs, but in general, all of the use cases involve spatial discretization and are typically based on variational formulations.

To simplify the notation, we will rewrite the reduced system as:

$$ g(u_c) + F_c(u_c,h(u_c)) = 0, $$

where $g(u_c)$ represents the time derivative. The precise form of h(𝑢𝑐) depends on the multiscale formulation, e.g., variational multiscale, homogenization, FE2, etc. In some cases, it may simply be a numerical approximation of $u_f$ given $u_c$, and in other cases it may be a model for the subgrid effects. However, in all cases we refer to h as the Macro-Micro-Macro map, since macro-scale information is sent to the local fine scale (subgrid) models, fine scale information is computed, and a reduced (upscaled) version of this information gets returned to the coarse scale problem (see the following figure).

MrHyDE is designed to utilize any type of subgrid model, but the primary demonstration is a multiscale Dirichlet-to-Neumann map which is described in below. MrHyDE also requires the gradient of $h$ with respect to $u_c$ to enable accurate Jacobians for implicit time integration. This is often complicated since $h$ will can depend on $u_c$ either explicitly or implicitly through $u_f$, which depends on $u_c$. Thus, the subgrid models typically compute

$$ \frac{\partial h}{\partial u_c} = \frac{\partial h}{\partial u_c} + \frac{\partial h}{\partial u_f}\frac{\partial u_f}{\partial u_c}. $$

Computing $\partial u_f / \partial u_c$ is the most challenging step. We first describe it for the steady-state case before moving on to the transient case.

Steady-state Case

In this section, we focus on solving the following steady-state multiscale system,

$$ F_c(u_c,u_f) = 0, $$ $$ F_f(u_c,u_f) = 0, $$

which is a simplification of the original coupled system. Recall that we assume that the fine scale equation can be localized, so given the current approximation for $u_c$, we can solve

$$F_f(u_c,u_f) = 0$$

for $u_f$. Note that this does assume that the fine scale solution can be uniquely determined given any coarse scale data, but this is a fairly reasonable assumption. The goal is then to solve the nonlinear coarse scale equation,

$$ F_c(u_c,h(u_c)) = 0, $$

for $u_c$. The default procedure in MrHyDE is to use Newton’s method, which, given an approximation $u_c^{(k)}$, solves

$$ J_c^{(k)} \Delta_c = -F_c\left(u_c^{(k)},h(u_c^{(k)})\right) $$

where

$$ J_c^{(k)} = \frac{\partial F_c}{\partial u_c}\big|_{u_c^{(k)}}, $$

and we set $_c^{(k+1)} = u_c^{(k)} +\Delta_c$. The iterations continue until a user-prescribed convergence 𝑐𝑐𝑐 tolerance, typically a relative tolerance on the norm of residual, is met. For the sake of simplicity, we will drop the iteration index on 𝑢𝑐 in the remainder of this section. In order to evaluate the residual, we need to be able to evaluate the Macro-Micro-Macro map, $h(u_c)$. This typically requires the following steps:

1. Either project the coarse scale data, $u_c$, onto the fine scale mesh or simply evaluate $u_c$ at the fine scale quadrature points.
2. Use another Newton-based procedure to solve for $u_f$ . Given an approximation $u_f^{(i)}$ (noting that $i$ is different from the outer Newton solve index $k$), we solve

$$ J_f^{(i)} \Delta_f = -F_f\left(u_f^{(i)},u_c)\right), $$

where

$$ J_f^{(i)} = \frac{\partial F_f}{\partial u_f}\big|_{u_f^{(i)}}, $$

and we set $u_f^{(i+1)} = u_f^{(i)} +\Delta_i$. Continue until convergence criteria is met.

3. Project $u_f$ back to the coarse scale or directly use the approximation.

At the discrete level, the projection steps above may involve a matrix-vector product or integrating the solutions against the appropriate test functions. MrHyDE precomputes the coarse scale basis functions at the fine scale quadrature points, so all of these integrals occur at the fine scale.

The coarse scale Jacobian can be written as,

$$ J_c = \frac{\partial}{\partial u_c}F_c(u_c,h(u_c)) = \frac{\partial F_c}{\partial u_c} + \frac{\partial F_c}{\partial h} \frac{\partial h}{\partial u_c}. $$

Computing $\partial F_c/\partial h$ is typically straightforward, but computing $\partial h/\partial u_c$, i.e., the derivative of the upscaled information with respect to the coarse scale data, requires more work. MrHyDE assumes that the subgrid models will provide $\partial h/\partial u_c$, which can be written as

$$ \frac{\partial h}{\partial u_c} = \frac{\partial h}{\partial u_f}\frac{\partial u_f}{\partial u_c}, $$

The key is computing $\partial u_f/\partial u_c$. To obtain this information, we differentiate the final scale equation with respect to $u_c$,

$$ \frac{\partial}{\partial u_c} F_f(u_c,u_f) = \frac{\partial F_f}{\partial u_c} + \frac{\partial F_f}{\partial u_f}\frac{\partial u_f}{\partial u_c} = 0. $$

We note that $\partial F_f/\partial u_f = J_f$, which we have already computed. In MrHyDE, we store the LU factorization (if using a direct method) or multilevel preconditioner (if using an iterative method) for the last Jacobian in the subgrid Newton solver, and solve

$$ J_f \frac{\partial u_f}{\partial u_c} = -\frac{\partial F_f}{\partial u_c}, $$

which typically involves one additional linear solve per degree-of-freedom in $u_c$. Since the matrix is already assembled and numerically factorized, this sensitivity propagation is typically faster than the actual subgrid nonlinear solve for $u_f$.

Transient Case

For transient simulations, we consider separately the cases where the nonlinear elimination occurs before or after the discretization in time. If we first discretize the coupled system in time, e.g., with a diagonally implicit RK scheme, and then perform nonlinear elimination to compute each of the stage solutions, then we call this synchronous time integration. Conversely, if the coarse scale and fine scale equations evolve at different time scales, then we might want to pursue a multirate formulation where the fine scale equation uses a different time integrator than the coarse scale equation. We refer to this as asynchronous time integration.

Synchronous Time Integration

First, we consider the simpler case where we discretize in time and then perform nonlinear elimination. Note that this does imply that both the coarse and fine scale equations evolve at the same rate. Recall that in MrHyDE, time integration schemes are implemented using a Butcher tableau and a BDF formula. Once discretized, we solve a sequence of nonlinear systems

$$ g(z_{c,i}) + F_c(m(z_{c,i}),m(z_{f,i})) = 0, $$ $$ g(z_{f,i}) + F_f(m(z_{c,i}),m(z_{f,i})) = 0, $$

for the stage solutions $z_{f,i}$ and $z_{c,i}$, where $g$ and $m$ involve the current, and possibly previous, stage solutions, and the step solutions are computed using

$$ u_{f,n+1} = u_{f,n} +\sum_{i=1}^s \left(z_{f,i} - u_{f,n} \right), $$ $$ u_{c,n+1} = u_{c,n} +\sum_{i=1}^s \left(z_{c,i} - u_{c,n} \right), $$

The reduced version of this problem solves

$$ g(z_{c,i}) + F_c(m(z_{c,i}),h(m(z_{c,i}))) = 0,$$

for $z_{c,i}$.

Asynchronous Time Integration

The asynchronous case is more complicated due to the fact that the coarse and fine scale integrators each have their own time step size, BDF formula and Butcher tableau. In addition, if the fine scale integrator uses a finer step size, then the coarse scale solution needs to be interpolated in time. For simplicity, we describe the solution process assuming we are using a Backward Euler integrator at both the coarse and fine scale, with $M$ substeps at the fine scale. Other variations are implemented and the solution process is similar.

Given $u_{c,n}$, the reduced version of the problem solves

$$ \frac{u_{c,n+1}-u_{c,n}}{\Delta t_c} + F_c(u_{c,n+1},h(u_{c,n+1})) = 0, $$

where we have used the fact that the backward Euler scheme is fairly simple to omit $g$ and $m$. Assuming that $u_{f,n}$ is given, we need to time step the fine scale model to obtain $u_{f,n+1}$ and $\partial u_{f,n+1}/\partial u_{c,n+1}$. We obtain this by taking $M$ subgrid time steps to go from $t_n$ to $t_{n+1}$. To simplify the notation, we set $u_f^{(0)} = u_{f,n}$ and compute

$$ \frac{u_{f}^{(1)} - u_{f}^{(0)}}{\Delta t_f} + F_f(p(t^{(1)},u_{c,n+1}),u_{f}^{(1)}) = 0 $$ $$ \frac{u_{f}^{(2)} - u_{f}^{(1)}}{\Delta t_f} + F_f(p(t^{(2)},u_{c,n+1}),u_{f}^{(2)}) = 0 $$ $$ \frac{u_{f}^{(3)} - u_{f}^{(2)}}{\Delta t_f} + F_f(p(t^{(3)},u_{c,n+1}),u_{f}^{(3)}) = 0 $$ $$ \frac{u_{f}^{(4)} - u_{f}^{(3)}}{\Delta t_f} + F_f(p(t^{(4)},u_{c,n+1}),u_{f}^{(4)}) = 0 $$ $$ \vdots \hspace{3cm} \quad \vdots $$

where $t^{(k)} = t^{(k)} + k\Delta t_f$ and $p(t^{(k)}, u_{c,n+1})$ denotes the interplant of $u_{c,n+1}$ in time at $t^{(k)}$. At the end, we set $u_{f,n+1} = u_f^{(M)}$. To compute $\partial u_{f,n+1} / \partial u_{c,n+1}$, we need to propagate this through the solution sequence appropriately. For ease of notation, set $v^{(k)} = \partial u_{f}^{(k)}/\partial u_{c,n+1}$. Then we have

$$ \frac{v^{(1)} - v^{(0)}}{\Delta t_f} + \frac{\partial F_f}{\partial p}\frac{\partial p}{\partial u_{c,n+1}} + \frac{\partial F_f}{\partial u_f}v^{(1)} = 0 $$ $$ \frac{v^{(2)} - v^{(1)}}{\Delta t_f} + \frac{\partial F_f}{\partial p}\frac{\partial p}{\partial u_{c,n+1}} + \frac{\partial F_f}{\partial u_f}v^{(2)} = 0 $$ $$ \frac{v^{(3)} - v^{(2)}}{\Delta t_f} + \frac{\partial F_f}{\partial p}\frac{\partial p}{\partial u_{c,n+1}} + \frac{\partial F_f}{\partial u_f}v^{(3)} = 0 $$ $$ \frac{v^{(4)} - v^{(3)}}{\Delta t_f} + \frac{\partial F_f}{\partial p}\frac{\partial p}{\partial u_{c,n+1}} + \frac{\partial F_f}{\partial u_f}v^{(4)} = 0 $$ $$ \vdots \hspace{3cm} \quad \vdots $$

where $v^{(0)} = 0$ and we set $\partial u_{f,n+1}/\partial u_{c,n+1} = v^{(M)}$.