Matrix-free Approach Research Articles

Investigation of biomolecular dynamics by sensitivity-enhanced 1H–2H CPMAS NMR using matrix-free dynamic nuclear polarization

Molecular dynamics of functional groups contain valuable information about structural properties and functional activities in biomolecules. NMR spectroscopy is a sensitive tool for the investigation of molecular dynamics over a wide range of timescales and thus may deepen the understanding of the biomolecules of interest. Here, we present an approach to use DNP-enhanced 2H NMR to study dynamics of selectively deuterated methyl groups in insoluble proteins such as amyloid beta (Aβ) fibrils. We adopted and optimized the matrix-free DNP approach by varying the amount of added polarizing agent as well as the rehydration level of model proteins. We show that the DNP enhancement obtained in 1H–2H cross-polarization (CP) MAS spectra may increase the sensitivity for selectively deuterated Aβ fibril samples by more than one order of magnitude, accelerating the collection of spin-lattice relaxation data in the DNP-accessible temperature range between 100 and 150 K by up to 400-fold. However, below the coalescence temperature, which describes the transition from the fast to the slow exchange regime, the experimentally obtained relaxation time constants suffer from a paramagnetic relaxation enhancement effect due to the presence of the polarizing agent. This seems to be a general effect for biomolecules as it is also confirmed for two other protein model systems. Our demonstration opens the possibility to extend the scope of 2H NMR for dynamics measurements to effective concentrations and/or temperatures below what is currently accessible; however, the observed interplay between paramagnetic relaxation and molecular dynamics also emphasizes the necessity for a better understanding of relaxation effects in DNP-enhanced NMR.

A matrix-free parallel two-level deflation preconditioner for two-dimensional heterogeneous Helmholtz problems

We propose a matrix-free parallel two-level deflation method combined with the Complex Shifted Laplacian Preconditioner (CSLP) for two-dimensional heterogeneous Helmholtz problems encountered in seismic exploration, antennas, and medical imaging. These problems pose challenges in terms of accuracy and convergence due to scalability issues with numerical solvers. Motivated by the limitations imposed by excessive computational time and memory constraints when employing a sequential solver with constructed matrices, we parallelize the two-level deflation method without constructing any matrices. Our approach utilizes preconditioned Krylov subspace methods and approximates the CSLP preconditioner with a parallel geometric multigrid V-cycle. For the two-level deflation, standard inter-grid deflation vectors and further high-order deflation vectors are considered. As another main contribution, the matrix-free Galerkin coarsening approach and a novel re-discretization scheme as well as high-order finite-difference schemes on the coarse grid are studied to obtain wavenumber-independent convergence. The optimal settings for an efficient coarse-grid problem solver are investigated. Numerical experiments of model problems show that the wavenumber independence has been obtained for medium wavenumbers. The matrix-free parallel framework shows satisfactory weak and strong parallel scalability.

A Coupling Physics Model for Real-Time 4D Simulation of Cardiac Electromechanics

Cardiac simulators can assist in the diagnosis of heart disease and enhance human understanding of this leading cause of mortality. The coupling of multiphysics, such as electrophysiology and active–passive mechanics, in the simulation of the heart poses challenges in utilizing existing methodologies for real-time applications. The low efficiency of physically-based simulation is mostly caused by the need for electrical-stress conduction to use tiny time steps in order to prevent numerical instability. Additionally, the mechanical simulation experiences sluggish convergence when dealing with significant deformation and stiffness, and there are also concerns regarding volume inversion. We provide a coupling physics model that transforms the active–passive dynamics into multiphysics solving constraints, aiming at boosting the real-time efficiency of the cardiac electromechanical simulation. The multiphysics processes are initially divided into two levels: cell-level electrical stimulation and organ-level electrical-stress diffusion/conduction. This separation is achieved by employing operator splitting in combination with the quasi-steady-state method, which simplifies the system equations. Next, utilizing spatial discretization, we employ the matrix-free conjugate gradient approach to solve the electromechanical model, therefore improving the efficiency of the simulation. The experimental results illustrate that our simulation model is capable of replicating intricate cardiac physiological phenomena, including 3D spiral waves and rhythmic contractions. Moreover, our model achieves a significant advancement in real-time computation while maintaining a comparable level of accuracy to current methods. This improvement is advantageous for interactive medical applications.

Tensor networks for solving the time-independent Boltzmann neutron transport equation

Tensor network techniques, known for their low-rank approximation ability that breaks the curse of dimensionality, are emerging as a foundation of new mathematical methods for ultra-fast numerical solutions of high-dimensional Partial Differential Equations (PDEs). Here, we present a mixed Tensor Train (TT)/Quantized Tensor Train (QTT) approach for the numerical solution of time-independent Boltzmann Neutron Transport equations (BNTEs) in Cartesian geometry. Discretizing a realistic three-dimensional (3D) BNTE by (i) diamond differencing in space, (ii) multigroup-in-energy, and (iii) discrete ordinates collocation in angle leads to large generalized eigenvalue problems that generally require a matrix-free approach and large computer clusters. Starting from this discretization, we construct a TT representation of the PDE fields and discrete operators, followed by a QTT representation of the TT cores. We then solve the tensorized generalized eigenvalue problem using a fixed-point scheme with tensor network optimization techniques. We validate our approach by applying the method to two examples of 3D neutron transport problems, currently solved by the Los Alamos National Laboratory PARallel TIme-dependent SN (PARTISN) solver.1 We demonstrate that our TT/QTT method, executed on a standard desktop computer, leads to large compression. This allows for the storage of terrabyte-sized neutron angular flux eigenvectors in megabytes. Additionally, we create megabyte-sized full access TT representations of yottabyte-sized transport matrix operators. By leveraging the TT operators and solution methods, we obtain a 7500 times speedup when compared to the PARTISN solution time with an error of less than 10−5.

Model-Based 3-D X-Ray Induced Acoustic Computerized Tomography.

X-ray-induced acoustic (XA) computerized tomography (XACT) is an evolving imaging technique that aims to reconstruct the X-ray energy deposition from XA measurements. Main challenges in XACT are the poor signal-to-noise ratio and limited field-of-view, which cause artifacts in the images. We demonstrate the efficacy of model-based (MB) algorithms for three-dimensional XACT and compare with the traditional algorithms. The MB algorithm is based on iterative, matrix-free approach for regularized-least-squares minimization corresponding to XACT. The matrix-free-LSQR (MF-LSQR) and the non-iterative model-backprojection (MBP) reconstructions were evaluated and compared with universal backprojection (UBP), time-reversal (TR) and fast-Fourier transform (FFT)-based reconstructions for numerical and experimental XACT datasets. The results demonstrate the capability of MF-LSQR algorithm to reduce noisy artifacts thus yielding better reconstructions. MBP and MF-LSQR algorithms perform particularly well with the experimental XACT dataset, where noise in signals significantly affects the reconstruction of the target in UBP and FFT-based reconstructions. The TR reconstruction for experimental XACT are comparable to MF-LSQR, but takes thrice as much time and filters the frequency components greater than maximum frequency supported by the grid, resulting loss of resolution. The MB algorithms are able to overcome the challenges in XACT and hence are vital for the clinical translation of XACT.

Krylov Methods for Large-Scale Dynamical Systems: Application in Fluid Dynamics

AbstractIn fluid dynamics, predicting and characterizing bifurcations, from the onset of unsteadiness to the transition to turbulence, is of critical importance for both academic and industrial applications. Different tools from dynamical systems theory can be used for this purpose. In this review, we present a concise theoretical and numerical framework focusing on practical aspects of the computation and stability analyses of steady and time-periodic solutions, with emphasis on high-dimensional systems such as those arising from the spatial discretization of the Navier–Stokes equations. Using a matrix-free approach based on Krylov methods, we extend the capabilities of the open-source high-performance spectral element-based time-stepper Nek5000. The numerical methods discussed are implemented in nekStab, an open-source and user-friendly add-on toolbox dedicated to the study of stability properties of flows in complex three-dimensional geometries. The performance and accuracy of the methods are illustrated and examined using standard benchmarks from the fluid mechanics literature. Thanks to its flexibility and domain-agnostic nature, the methodology presented in this work can be applied to develop similar toolboxes for other solvers, most importantly outside the field of fluid mechanics.

Scalable matrix-free solver for 3D transfer of polarized radiation in stellar atmospheres

We present an efficient and massively parallel solution strategy for the transfer problem of polarized radiation, for a 3D stationary medium out of local thermodynamic equilibrium. Scattering processes are included accounting for partial frequency redistribution effects. Such a setting is one of the most challenging ones in radiative transfer modeling. The problem is formulated for a two-level atomic model, which allows linearization. The discrete ordinate method alongside an exponential integrator are used for discretization. Efficient solution is obtained with a Krylov method equipped with a tailored physics-based preconditioner. A matrix-free approach results in a lightweight implementation, suited for tackling large problems. Near-optimal strong and weak scalability are obtained with two complementary decompositions of the computational domain. The presented approach made it possible to perform simulations for the Ca i line at 4227Å with more than 109 degrees of freedom in less than half an hour on massively parallel machines, always converging in a few iterations for the proposed tests.

Limits of entrainment of circadian neuronal networks.

Circadian rhythmicity lies at the center of various important physiological and behavioral processes in mammals, such as sleep, metabolism, homeostasis, mood changes, and more. Misalignment of intrinsic neuronal oscillations with the external day-night cycle can disrupt such processes and lead to numerous disorders. In this work, we computationally determine the limits of circadian synchronization to external light signals of different frequency, duty cycle, and simulated amplitude. Instead of modeling circadian dynamics with generic oscillator models (e.g., Kuramoto-type), we use a detailed computational neuroscience model, which integrates biomolecular dynamics, neuronal electrophysiology, and network effects. This allows us to investigate the effect of small drug molecules, such as Longdaysin, and connect our results with experimental findings. To combat the high dimensionality of such a detailed model, we employ a matrix-free approach, while our entire algorithmic pipeline enables numerical continuation and construction of bifurcation diagrams using only direct simulation. We, thus, computationally explore the effect of heterogeneity in the circadian neuronal network, as well as the effect of the corrective therapeutic intervention of Longdaysin. Last, we employ unsupervised learning to construct a data-driven embedding space for representing neuronal heterogeneity.

Magnetic Nanoparticle-Based Surface-Assisted Laser Desorption/Ionization Mass Spectrometry for Cosmetics Detection in Contaminated Fingermarks: Magnetic Recovery and Surface Roughness.

In this work, we propose a matrix-free approach for the analysis of fingermarks (FMs) contaminated with five cosmetic products containing different active pharmaceutical ingredients (APIs) using surface-assisted laser desorption/ionization mass spectrometry (SALDI-MS). For this purpose, a magnetic SALDI substrate based on Fe3O4-CeO2 magnetic nanoparticles was prepared, characterized, and optimized for the analysis of contaminated FMs without sample pretreatment. Initially, groomed FM and cosmetic products were separately analyzed, and their major components were successfully detected. Subsequently, FMs contaminated with Ordinary serum and Skinoren, Dermovate, Bepanthen, and Eucerin creams were analyzed, and components of FM and cosmetics were detected. The stability of the cosmetics in FMs was studied over an interval of 28 days, and all components showed good stability in FM for 4 weeks. Recovery of contaminated FMs from different surfaces utilizing a few microliters of the magnetic substrate was carried out using a simple external magnetic field from ceramic, plastic, metal, and glass. Successful retrieval of the API and FM components was achieved with magnetic recovery, and glass exhibited the best recovery, whereas ceramic tile demonstrated the lowest recovery. This was supported by atomic force microscopy study, which revealed that the ceramic surface had higher roughness than the other surfaces employed in this study, which adversely affected the magnetic maneuvering. This proof-of-concept investigation extends the application of SALDI-MS in forensic analysis of contaminated FMs by exploring cosmetics as exogenous materials and their stability and recovery from different surfaces.

Implicit two-derivative deferred correction time discretization for the discontinuous Galerkin method

In this paper, we use an implicit two-derivative deferred correction time discretization approach and combine it with a spatial discretization of the discontinuous Galerkin spectral element method to solve (non-)linear PDEs. The resulting numerical method is high order accurate in space and time. As the novel scheme handles two time derivatives, the spatial operator for both derivatives has to be defined. This results in an extended system matrix of the scheme. We analyze this matrix regarding possible simplifications and an efficient way to solve the arising (non-)linear system of equations. It is shown how a carefully designed preconditioner and a matrix-free approach allow for an efficient implementation and application of the novel scheme. For both, linear advection and the compressible Euler equations, up to eighth order of accuracy in time is shown. Finally, it is illustrated how the method can be used to approximate solutions to the compressible Navier-Stokes equations.

Matrix-free approaches for GPU acceleration of a high-order finite element hydrodynamics application using MFEM, Umpire, and RAJA

With the introduction of advanced heterogeneous computing architectures based on GPU accelerators, large-scale production codes have had to rethink their numerical algorithms and incorporate new programming models and memory management strategies in order to run efficiently on the latest supercomputers. In this work we discuss our co-design strategy to address these challenges and achieve performance and portability with MARBL, a next-generation multi-physics code in development at Lawrence Livermore National Laboratory. We present a two-fold approach, wherein new hardware is used to motivate both new algorithms and new abstraction layers, resulting in a single source application code suitable for a variety of platforms. Focusing on MARBL’s ALE hydrodynamics package, we demonstrate scalability on different platforms and highlight that many of our innovations have been contributed back to open-source software libraries, such as MFEM (finite element algorithms) and RAJA (kernel abstractions).

Enhancing scalability of a matrix-free eigensolver for studying many-body localization

We propose several techniques to enhance the parallel scalability of a matrix-free eigensolver designed for studying many-body localization (MBL) of quantum spin chain models with nearest-neighbor interactions and on-site disorder. This type of problem is computationally challenging because the dimension of the associated Hamiltonian matrix grows exponentially with respect to the number of spins L, and we need to average over different realizations of the random disorder to obtain relevant statistical behavior. For each disorder realization, we need to compute eigenvalues from different regions of the spectrum and their corresponding eigenvectors. In previous work, the interior eigenstates for a single eigenvalue problem are computed via the shift-and-invert Lanczos algorithm. Due to the extremely high memory footprint of the LU factorizations, this technique is not well suited for large L’s. For example, we need thousands of compute nodes on modern high performance computing infrastructures to go beyond L = 24. The matrix-free approach does not suffer from this memory bottleneck, however, its scalability is limited by a computation and communication load imbalance. To reduce this imbalance and to significantly enhance the scalability of the matrix-free eigensolver, we reorder the matrix and leverage the consistent space runtime, CSPACER. We also show its efficiency in managing irregular communication patterns at scale compared to optimized MPI non-blocking two-sided and one-sided RMA implementation variants. This effort enables us to study MBL for spin chains with a larger number of spins. The efficiency and effectiveness of the proposed algorithm is demonstrated by computing eigenstates on a massively parallel many-core high performance computer.

A fast matrix-free approach to the high-order control volume finite element method with application to low-Mach flow

A fast matrix-free formulation of the control volume finite element method is presented, requiring much less memory and computational work than previous efforts. The method is implemented and evaluated as a solver for low-Mach flow, including the evaluation of a preconditioning strategy for the pressure Poisson equation. The efficiency and scaling with polynomial order is evaluated on simple turbulent flows of interest, with appropriate solution quality metrics, and compared with a reference node-centered finite volume discretization. For a turbulent channel flow test, we show improvement in computational work for a given accuracy with the high-order scheme. The performance on a GPU accelerated platform is also investigated, with benefit shown for the matrix-free discretization.

Tensorized low-rank circulant preconditioners for multilevel Toeplitz linear systems from high-dimensional fractional Riesz equations

The solution of d-dimensional fractional partial differential equations (FPDEs) is challenging in numerical computation. This paper proposes a low-rank preconditioned conjugate gradient (PCG) method for discretizations with a low-rank right-hand side stemming from the d-dimensional fractional Riesz equation. At each iteration step, a combination of low-rank in Tensor Train (TT) format technique and Fast Fourier Transform (FFT) is used to accelerate multilevel Toeplitz matrix-vector multiplications. Furthermore, a low-rank circulant preconditioner is constructed based on the spectrum of the coefficient matrix. The reported numerical results demonstrate the competitiveness of the new method compared with a state-of-the-art matrix-free approach based on FFT.

Development of efficient SALDI substrate based on Au–TiO2 nanohybrids for environmental and forensic detection of dyes and NSAIDs

Herein, a matrix-free approach is presented for comprehensive environmental and forensic analysis of dyes and nonsteroidal anti-inflammatory drugs (NSAIDs) using Au–TiO2 nanohybrids coupled with surface-assisted pulsed laser desorption ionization-mass spectrometry (SALDI-MS). The Au–TiO2 nanohybrids was prepared and characterized using inductively coupled plasma-optical emission spectrometry (ICP-OES), X-ray photoelectron spectroscopy (XPS), X-ray diffraction (XRD), surface area measurements, ultraviolet–visible (UV–vis) spectroscopy, transmission electron microscopy (TEM), and scanning electron microscopy-energy dispersive spectrometry (SEM-EDS). Initially, the optimal Au content was assessed using the survival yield (SY) method, confirming that 7.5% Au content on the TiO2 surface offered the highest ionization efficiency. Subsequently, environmental analyses of dyes and NSAIDs in water samples were performed, and sensitive detection of all analytes was achieved with limits of detection (LODs) ranging from 10.0 ng mL−1 to 10.0 fg mL−1 and good spot-to-spot reproducibility. Additionally, the effect of potential contaminants commonly found in environmental samples, such as salts, surfactants and pesticides was also considered. Despite signal intensity reduction at high concentrations of some salts, the target analytes were detected, while the presence of surfactants and pesticides did not cause significant signal intensity reduction. Additionally, dyed and undyed Tetoron fibers and the effect of adhesive tape were evaluated. Direct analysis of the dyed Tetoron fibers on the target plate, using the nanohybrids, enabled higher detection sensitivity of the dyes, in addition to adducts of polystyrene and cellulose, the main components of the fiber. Finally, NSAIDs in oral fluid were analyzed and sensitive detection of the analytes was observed using the nanohybrids with LODs and LOQs in the range of 0.1–10 ng mL−1 and 1–20 ng mL−1, respectively. The trueness of the exact mass was in the range of 0.64–6.2 ppm while the recovery of the spiked samples was in the range of 82.90–107.54%% indicating the efficiency of the Au–TiO2 nanohybrids as SALDI substrate. Thus, the Au–TiO2 nanohybrids hold considerable promise in terms of sensitivity, reproducibility, and LOD, and may significantly contribute to environmental and forensic identification.

Matrix-free TriGlobal adjoint stability analysis of compressible Navier–Stokes equations

A numerical method for TriGlobal (i.e., fully three-dimensional) adjoint stability analysis for compressible flows was developed and is presented in this paper. The developed method solves the adjoint stability problem using a matrix-free method based on Krylov-Schur method and a time-stepping approach. Because of the low memory (RAM) requirement of the matrix-free approach, the developed method can analyze fully three-dimensional flows that are difficult to analyze with conventional matrix-based methods. To perform time-stepping on the adjoint variables, the adjoint equations, including appropriate boundary conditions for the compressible Navier–Stokes equations, were derived. The equations were discretized using a finite compact difference method, and time integration was conducted using the three-step third-order Runge–Kutta method. A flow field around a two-dimensional square cylinder and a cubic cavity flow were analyzed using the developed method, and it was confirmed that the method reproduces the dominant instabilities reported in the literature. In addition, in the square cylinder flow analysis, the receptivity and sensitivity regions of the secondary wake mode, which corresponds to far wake instability, were clarified. Finally, the TriGlobal direct and adjoint stabilities of compressible flows over a finite width cavity were analyzed for the first time, and it was proven that large-scale adjoint stability analysis can be performed with the developed method. The results also show that instability phenomena similar to those obtained with BiGlobal stability analysis appear, but sidewall effects exist.

Induced polarization of volcanic rocks. 4. Large-scale induced polarization imaging

SUMMARY Thanks to the emergence of new technologies developed with the goal of performing large-scale galvanometric induced polarization surveys and thanks a better understanding of the underlying physics of induced polarization, this geophysical method can now be applied in the field of volcanology and geothermal resources assessment. A new approach is developed here for directly inverting the primary and secondary electric fields recorded at a set of independent stations when injecting a primary current. The use of independent stations to measure the primary and secondary electrical fields improves the quality of the data by reducing the capacitive coupling effects inherent to systems based on long cables. It avoids issues associated with using the same electrodes for both current injection and voltage measurements and negative apparent resistivity and chargeability values. With such acquisitions, we can perform true 3-D surveys in areas characterized by complex topography such as volcanoes. The numerical scheme we developed returns as output the electrical conductivity and chargeability fields. The implemented methodology presents several advantages. The first is the use of data types at the stations, for example the electric field intensity, that are independent from the local geometrical station parameters such as electrode spacing and dipole orientation. The second advantage lies in the suitability of the proposed approach to perform large-scale applications since we use a matrix-free approach that does not require the assembly of the Jacobian matrices. The third concerns the possibility of performing the inversion on complex geometries through a consistent use of the finite element method on unstructured meshes in combination with algebraic multigrid preconditioning for the regularization and the solution of the forward and adjoint problems. The computation of 3-D sensitivity maps can also be a real asset in survey design. After validating our approach with a benchmark synthetic case study, we test it on a large-scale induced polarization survey that mimic true field conditions on a volcanic environment with rough topography. Our tests demonstrate the high potential of this electric field approach in volcanology especially for deep (3 km) imagining of the internal structure of volcanoes, which in turn could improve our understanding of hydrothermal systems and allow the monitoring of active volcanoes and the potential risk of collapse.

Matrix-free multigrid solvers for phase-field fracture problems

In this work, we present a framework for the matrix-free solution to a monolithic quasi-static phase-field fracture model with geometric multigrid methods. Using a standard matrix-based approach within the Finite Element Method requires lots of memory, which eventually becomes a serious bottleneck. A matrix-free approach overcomes this problem and greatly reduces the amount of required memory, allowing to solve larger problems on available hardware. One key challenge is concerned with the crack irreversibility for which a primal–dual active set method is employed. Here, the active set values of fine meshes must be available on coarser levels of the multigrid algorithm. The developed multigrid method provides a preconditioner for a generalized minimal residual (GMRES) solver. This method is used for solving the linear equations inside Newton’s method for treating the overall nonlinear-monolithic discrete displacement/phase-field formulation. Several numerical examples demonstrate the performance and robustness of our solution technology. Mesh refinement studies, variations in the phase-field regularization parameter, iterations numbers of the linear and nonlinear solvers, and some parallel performances are conducted to substantiate the efficiency of the proposed solver for single fractures, multiple pressurized fractures, and a L-shaped panel test in three dimensions.

A scalable matrix-free spectral element approach for unsteady PDE constrained optimization using PETSc/TAO

We provide a new approach for the efficient matrix-free application of the transpose of the Jacobian for the spectral element method for the adjoint-based solution of partial differential equation (PDE) constrained optimization. This results in optimizations of nonlinear PDEs using explicit integrators where the integration of the adjoint problem is not more expensive than the forward simulation. Solving PDE constrained optimization problems entails combining expertise from multiple areas, including simulation, computation of derivatives, and optimization. The Portable, Extensible Toolkit for Scientific computation (PETSc) together with its companion package, the Toolkit for Advanced Optimization (TAO), is an integrated numerical software library that contains an algorithmic/software stack for solving linear systems, nonlinear systems, ordinary differential equations, differential algebraic equations, and large-scale optimization problems and, as such, is an ideal tool for performing PDE-constrained optimization. This paper describes an efficient approach in which the software stack provided by PETSc/TAO can be used for large-scale nonlinear time-dependent problems. Time integration can involve a range of high-order methods, both implicit and explicit. The PDE-constrained optimization algorithm used is gradient-based and seamlessly integrated with the simulation of the physical problem.

Adjoint-based inversion for porosity in shallow reservoirs using pseudo-transient solvers for non-linear hydro-mechanical processes

Porous flow is of major importance in the shallow subsurface, since it directly impacts on reservoir-scale processes such as waste fluid sequestration or oil and gas exploration. Coupled and non-linear hydro-mechanical processes describe the motion of a low-viscous fluid interacting with a higher viscous porous rock matrix. This two-phase flow may trigger the initiation of solitary waves of porosity, further developing into vertical high-porosity pipes or chimneys. These preferred fluid escape features may lead to localised and fast vertical flow pathways potentially problematic in the case of for instance CO2 sequestration. Constraining the porosity and the non-linearly related permeability distribution in such environments is a major challenge. Although seismic imaging methods accurately localise the high-porosity chimneys in the inverted wave-speed field, the conversion to porosity is not straightforward. We develop an inversion framework to reconstruct the unknown porosity field using relevant observable quantities such as subsurface fluid fluxes. We introduce the adjoint framework for the two-phase flow equations, which allows for efficient calculations of the pointwise gradients of the flow solution with respect to the porosity. We then use the gradients in a gradient descent method to invert for the pointwise porosity. We solve the forward and the adjoint equations using an iterative matrix-free pseudo-transient approach with the finite difference method. The proposed parallel solving procedure executes optimally on the latest many-core hardware accelerators such as GPUs. Numerical results show that an inversion for porosity is challenging if data is sparse since the porosity is very locally sensitive to the fluid flux. We introduce the concepts of sensitivity kernels as employed in seismology for the set of two-phase equations and suggest this approach as a standard for future studies.