-  -   home   -  -  -   research   -  -  -   students   -  -  -   courses   -  -  -   publications   -  -  -   vitae   -  -  

Research


Flux partitioning for air-sea interaction
Computationally demanding applications like climate prediction and weather forecasting often
require air and sea codes to be coupled together. There are many complications. The air
and sea systems are quite different, so that their efficient simulation requires a heterogeneous
modeling approach. Various numerical methods are employed together simultaneously; in particular,
the air and sea model may use different spatial grid sizes and internal time steps. In order to
run these codes together, a popular approach is to use a so-called ''flux coupler'' that takes
in state and grid information from the two fluid modules and computes net fluxes of conserved
physical quantities. The fluxes are used to prescribe boundary conditions for the two fluid
codes. Among many issues that the flux coupler must handle, mismatching grids and time steps
introduce questions for constructing accurate, conservative and stable flux computations.

Over many years now, numerous papers have appeared that address various aspects of the time
stepping problem. Typical methods have a (formally) low order of accuracy in terms of error scaling with
respect to time step sizes. But methods development is nontrivial. Nonlinear fluxes and the related
boundary conditions have presented a formidable challenge to construct and analyze methods, leading to
some open theoretical problems. Conservation is important for fluxes, but numerous analyses indicate
a fundamental problem to construct methods that are both very stable and conservative. Connors, Howell
and Layton (SINUM, 2012) presented some non-conservative, low-order methods with one being unconditionally
stable. The stability of the other methods was left open and only some inconclusive conditional analysis
has appeared since for that issue. Based on these methods, it was shown for coupled fluids with
natural heat convection (Connors and Ganis, 2011) that numerical bias introduced through coupling errors
could potentially not be captured in parameter sensitivity studies, even when it outweighs other sources
of model uncertainty. The only known high-order method with unconditional stability
appeared recently (Aggul, Connors, Erkmen and Labovsky, SINUM, 2018). This approach uses the
deferred correction technique to iteratively achieve high-order accuracy from the stable, low-order
method of Connors, Howell and Layton, so it is still not conservative. The existence of a non-iterative,
unconditionally stable and high-order method also remains open, conservative or otherwise.

Also, the high computational cost for applications requires that the highly-optimized codes be used with the
flux-coupler connectivity, which imposes some constraints on the practical side for methods development. In
turn, there has been a need to develop theoretical coupling frameworks for a non-intrusive class of methods.
In lieu of developing methods with built-in unconditional stability properties, another approach has been to
investigate ways to mitigate stability issues. Some work has appeared exploring iteration between air and
sea modules over a time interval to tighten the coupling and improve accuracy conservatively (for example
see Lemarié, Blayo and Debreu, Procedia Computer Science, 2015). Many related coupling
configurations have recently been explored (Connors and Dolan, AAMM 2019) using least-squares based
flux computations to achieve high-order accuracy with moments of the flux conserved over a coupling time
interval. Experiments therein indicate that iteration can increase stability and accuracy, but not much
compared to decreasing time step sizes and using higher-order methods. Meanwhile, a natural question is
how to develop an adaptive method. A stability indicator that could be easily computed on-the-fly was
used by Connors and Dolan as a way to selectively iterate for stability, but some evidence shows that
adaptive time step selection is probably more useful. At present, adaptive methods for flux computations
are still in their infancy.

There have been further developments for air and sea codes as well as coupling software in recent years. As
next-generation computing platforms loom in the future, current code structures may become more adaptable and
computational details of model layering, mesh and time step parameters will likely change, in turn changing
the picture for coupling techniques. Methods that seem less applicable now may become more attractive later.


Numerical error estimation
Complex simulations have numerous sources of uncertainty that can make it difficult to interpret model
outputs. One source of uncertainty is numerical discretization error, which motivates techniques to compute
estimates of this error, in particular for use with computationally demanding simulations. This is a
different perspective from the classical techniques of numerical error estimation motivated primarily by
adaptive mesh refinement. One promising approach is the technique of nonlinear error transport, which
has had success even for computing errors near discontinuities, including shocks in fluids (see, for example,
Banks, Hittinger, Connors and Woodward; CMAME 2012 and Contemp. Math. 2013). This technique approximates
the entire field of errors for computed states in space and time, yielding quite detailed error information. It
is based on the idea of deriving the error equation corresponding to your governing equation, then modifying
it in a standardized way to make it amenable to subsequent numerical approximation, yielding an error estimate.

Another interesting approach is adjoint error estimation, which focuses instead on a functional
measurement of the error. An integral expression for the error in the functional measurement can be derived
by defining a certain adjoint problem, which depends on the primal (governing) problem and the choice of
functional. The integral expression is approximated by inserting computed states for the primal and adjoint
problems, yielding a computable error estimate. This approach has the advantage of being able to simultaneously
compute the error and adapt the computational grid according to the functional of interest. The accuracy of
the error estimate can again be quite good even with a discontinuous solution (Connors, Banks, Hittinger and
Woodward, SINUM 2013).

Interestingly, adjoint methods have also proved useful for separating discretization errors into
components, which is of particular use for multiscale and multiphysics applications. For example, in
multiphysics applications operator-splitting techniques are popular (Lie-Trotter, Strang, etc.) as a
way to resolve two (or more) different physical processes that happen simultaneously, by numerically separating
them to evolve independently for very short time intervals. In this way, the system can be simulated by
putting together individual codes that are highly optimized for the separate physical processes. However,
the artificial separation of the physical processes introduces a type of numerical error, which is coupled
to the underlying discretization errors from the code modules being used for each process. It raises a
question of the size of the so-called splitting error relative to the discretization errors. In fact,
separating errors into components and computing estimates of each component is often possible with sufficient
accuracy to distinguish their relative sizes well. Connors, Banks, Hittinger and Woodward demonstrated this
when splitting advection and diffusion operators (CMAME 2014), for example. Other sorts of error decomposition
is a subject of on-going research.

Computing dissipative solitons
In nature, there are many examples of pattern formation. Many such phenomena are modeled using reaction-diffusion
equations, in particular activator-inhibitor systems that allow a delicate balance between gain and loss of energy and
chemical concentrations. Against an ambient background state, localized structures like wave fronts and pulses are
commonly observed. Stable traveling pulses are often called dissipative solitons in this context, as they interact
with one another in particle-like manners. There is a voluminous body of literature both for mathematical analysis of
these patterns and numerical simulations to explore many of their properties.

Theoretical analysis has its limitations and many questions remain. In some particular parameter regimes and with
some systems, it is not known what sorts of patterns like standing or traveling waves may exist. It can be helpful
to explore by using numerical simulation, for which purpose the popular technique is parameter continuation. It may
be possible to calculate an initial wave profile in various ways for a particular choice of model parameters, then
incrementally adjust the parameters and update the computed profile to explore solutions along bifurcation curves
in parameter space. As an illustration, when seeking traveling dissipative solitons for certain systems, a
common approach has been to start from a standing wave solution (so wave speed zero). Continuation is then done
by slowly increasing the wave speed to find traveling solitons. For some problems, the result finds only
unstable waves, leading researchers to conjecture that traveling soliton solutions do not exist.

Recently, Choi and Connors proposed using a steepest-decent approach to compute solitons (see Choi and Connors, 2020).
It is a robust supplement for the classical continuation approach that applies when the governing equations have
a variational structure. In this case, the governing equations are the Euler-Lagrange equation for a certain variational
functional; minimizing this energy functional identifies a wave pattern. The method does not require any detailed
information about the solution a priori, so that one can search virtually anywhere in parameter space for solutions
independently, rather than approaching along a bifurcation curve (which may not connect to some solutions). In fact,
Choi and Connors found examples of dissipative solitons in 2D for a situation where they were believed not to exist. In
future studies, it will be interesting to see how the steepest-decent and continuation methods can be used together and
learn new things about patterns in nature.