Numerical experiments to assess robustness of variational data assimilation for hemodynamics

Research on mathematical modeling of hemodynamics has come a long way, with patient specific simulations based on vessel geometries gathered from CT or MR imaging now being fairly common in scientific journals. Besides geometry description, the flow rates and division of flow between vascular branches are arguably the most important patient specific quantity for quantification of hemodynamic forces. In pathological cases, flow rates may vary greatly over time and between patients and simulations based on reference values are therefore less realistic than they could be if flow rates were measured and included properly. Variational data assimilation (DA) is a promising technique for merging MRI velocity observations into simulations with the Navier-Stokes (NS) equations, resulting in a more trustworthy physical flow field for further analysis. DA for NS has received attention for decades in the context of geophysics and meteorology, however its application in hemodynamics applications is fairly recent (see Marta D'Elia, 2011). While the computational methods from previous research can be applied to a large degree, the problem specific characteristics in e.g. geophysics and hemodynamics differ in many respects. Good accuracy depends mainly on capturing the flow rates to a good accuracy, since flow in hemodynamics is otherwise largely determined by the pipe geometry. Hence it seems that DA is the best way to obtain good initial and boundary conditions for simulations. A good match between observed and simulated flow field will further serve to validate the modelling. The main challenges to overcome are accuracy and efficiency. Accuracy of DA is largely affected by the technological limits of imaging methods today, i.e. low resolution and noise, but also by the choice of computational method, which may involve a tradeoff between accuracy and efficiency. It is important to understand the characteristics of each method before we can perform larger studies and trust the results. The excessive computational cost for DA simulations with high resolution in time and space is a real limiting factor for research today, but much can be done on smaller models, lower resolution models, and in 2D before attacking the large scale problems. In particular, the study of methods in the light of synthetic data can give us a better understanding of how these models and methods behave in situations close to clinical data. In this study, we formulate the variational data assimilation problem as a constrained optimization problem, where the distance between the computed flow field and observed blood velocities is minimized, constrained by fluid flow physics goverened by the Navier-Stokes equations. The control parameters are the initial condition and boundary conditions at inlets and outlets, which are unknown in actual applications. The Navier-Stokes equations are discretized in space with the finite element method, and time-stepped using a splitting scheme to compute the transient solution efficiently. The L-BFGS optimization algorithm is employed, an algorithm which requires the computation of the gradient of the cost functional to minimize w.r.t. transient control functions. This gradient computation is achieved by solution of the adjoint equations of the time dependent problem. The finite element model of the state equations are implemented with FEniCS and the adjoint computation is automated by dolfin-adjoint. Three key enabling techniques for this automation are algorithmic differentiation, automated construction of the discrete adjoint equations, and optimal checkpointing. FEniCS contains a domain specific language with embedded algorithmic differentiation called UFL, which automates differentiation of the cost functional w.r.t. state and control variables, allowing the cost functional to be specified only once. Dolfin-adjoint records the forward model equation solves and automates replay and adjoint replay of the recorded equations. For large meshes and a large number of timesteps, not all intermediate solutions can be stored in memory or on disk, and an optimal checkpointing scheme from Revolve handles storage of a subset of the transient state vectors and subsequent recomputation of the forward problem during the adjoint computation. Our research goal is patient specific 3D simulations of transient blood flow in aneurysms, however we first need to develop stable methods and understand how the optimization methods, Navier-Stokes schemes, and physics of the optimization problem interact. In this talk we will present experiments with varying noise levels applied to artificial analytical cases with a range of Reynolds numbers, to assess the stability, accuracy and efficiency of the methods and problem setup. Steps will be taken towards more physiologically relevant test cases