Highly efficient NURBS-based isogeometric analysis for coupled nonlinear diffusion-reaction equations with and without advection

TL;DR

This work develops a high-order isogeometric analysis framework for coupled advection-diffusion-reaction systems by coupling semi-Lagrangian treatment of advection with a Strang operator split diffusion-reaction step, all discretized via NURBS-based IgA. The method avoids linearization of nonlinear reactions, uses a second-order BDF2 scheme for diffusion with $L^{2}$ projection to handle non-interpolatory NURBS, and employs explicit RK4 for reactions, yielding second-order accuracy and efficiency on complex geometries. Numerical experiments across fully nonlinear, exact-solution, Schnakenberg-Turing, and Gray-Scott problems demonstrate accurate pattern formation, geometry-driven effects on Turing patterns, and robustness to advection, highlighting the approach's potential for complex multiphysics in biological morphogenesis. The combination of $k$-refinement in IgA and operator splitting enables high accuracy with relatively coarse meshes, enabling efficient simulations on intricate domains. The results confirm the method's capability to capture pattern formation and the influence of geometry and transport on emergent phenomena.

Abstract

Nonlinear diffusion-reaction systems model a multitude of physical phenomena. A common situation is biological development modeling where such systems have been widely used to study spatiotemporal phenomena in cell biology. Systems of coupled diffusion-reaction equations are usually subject to some complicated features directly related to their multiphysics nature. Moreover, the presence of advection is source of numerical instabilities, in general, and adds another challenge to these systems. In this study, we propose a NURBS-based isogeometric analysis (IgA) combined with a second-order Strang operator splitting to deal with the multiphysics nature of the problem. The advection part is treated in a semi-Lagrangian framework and the resulting diffusion-reaction equations are then solved using an efficient time-stepping algorithm based on operator splitting. The accuracy of the method is studied by means of a advection-diffusion-reaction system with analytical solution. To further examine the performance of the new method on complex geometries, the well-known Schnakenberg-Turing problem is considered with and without advection. Finally, a Gray-Scott system on a circular domain is also presented. The results obtained demonstrate the efficiency of our new algorithm to accurately reproduce the solution in the presence of complex patterns on complex geometries. Moreover, the new method clarifies the effect of geometry on Turing patterns.

Figures (15)

• Figure 1: Mappings from the parent space through the parametric space to the physical space. For the purpose of illustration, we took the knot vectors $\Xi^{1} = \Xi^{2} = \bigl\{0,0,0,1/2,1,1,1\bigr\}$. Thus, the parametric space consists of four elements. We chose the degrees $p = q = 2$ which results in $16=4\times 4$ NURBS functions obtained as tensor product of the corresponding one-dimensional NURBS functions. The dashed blue line in the physical space refers to the control polygon associated to the circular object $\Omega$.
• Figure 2: A schematic diagram showing the main quantities used in the calculation of the departure points. Since NURBS functions lie in the parametric space, each quadrature point $(\bar{\xi_{g}},\bar{\eta_{g}})$ is first mapped from the parent space to the parametric space according to the mapping $\phi_{\mathcal{P}_{k}}$. This results in the point $\bm\xi_{k,g}=(\xi_{k,g},\eta_{k,g})$. The corresponding departure point $\mathcal{\widetilde{Y}}^{n}(\bm\xi_{k,g})$ is then calculated at the host element, see the dark element $\widetilde{\mathcal{P}}^{*}_{k}$, before being mapped to the physical element $\mathcal{P}_{k}$. This results in the departure point $\mathcal{Y}^{n}({\boldsymbol{x}}_{k,g})$. In this illustration, the local control points for NURBS functions (of degree 1) are indicated by small squares.
• Figure 3: Cross sections of the solution in Example \ref{['test0']} along the diagonal $x = y$, at time $t=1$ using different NURBS degrees, for two different mesh resolutions. Left: $32\times32$ elements. Right: $64\times64$ elements.
• Figure 4: Snap shots of the solution of Example \ref{['test1']} using different NURBS degrees. From left to right: $t = 0.5$, $t = 1$ and $t = 2$. The mesh resolution is fixed to $32\times32$ elements.
• Figure 5: Continuation of Fig. \ref{['Ex1_Fig1']}.