Engineering Papers⌕ Search

SEARCH · Engineering Papers

Results for “eigenvalue problems”

Search indexed NASA NTRS and DOE OSTI research on propulsion, heat transfer, battery materials and energy systems. Follow report and document links to the original sources.

Quote a phrase for an exact phrase match. Source license links do not imply unrestricted reuse.

At least 433 records · Page 24

Solving fractional Schrödinger-type spectral problems: Cauchy oscillator and Cauchy well

This paper is a direct offspring of the work of Garbaczewski and Stephanovich [“Lévy flights and nonlocal quantum dynamics,” J. Math. Phys. 54, 072103 (2013)] where basic tenets of the nonlocally induced random and quantum dynamics were analyzed. A number of mentions were made with respect to various inconsistencies and faulty statements omnipresent in the literature devoted to so-called fractional quantum mechanics spectral problems. Presently, we give a decisive computer-assisted proof, for an exemplary finite and ultimately infinite Cauchy well problem, that spectral solutions proposed so far were plainly wrong. As a constructive input, we provide an explicit spectral solution of the finite Cauchy well. The infinite well emerges as a limiting case in a sequence of deepening finite wells. The employed numerical methodology (algorithm based on the Strang splitting method) has been tested for an exemplary Cauchy oscillator problem, whose analytic solution is available. An impact of the inherent spatial nonlocality of motion generators upon computer-assisted outcomes (potentially defective, in view of various cutoffs), i.e., detailed eigenvalues and shapes of eigenfunctions, has been analyzed.

71 CLASSICAL AND QUANTUM MECHANICS, GENERAL PHYSIC↗

GW and Bethe-Salpeter study of small water clusters

We study within the GW and Bethe-Salpeter many-body perturbation theories the electronic and optical properties of small (H{sub 2}O){sub n} water clusters (n = 1-6). Comparison with high-level CCSD(T) Coupled-Cluster at the Single Double (Triple) levels and ADC(3) Green’s function third order algebraic diagrammatic construction calculations indicates that the standard non-self-consistent G{sub 0}W{sub 0}@PBE or G{sub 0}W{sub 0}@PBE0 approaches significantly underestimate the ionization energy by about 1.1 eV and 0.5 eV, respectively. Consequently, the related Bethe-Salpeter lowest optical excitations are found to be located much too low in energy when building transitions from a non-self-consistent G{sub 0}W{sub 0} description of the quasiparticle spectrum. Simple self-consistent schemes, with update of the eigenvalues only, are shown to provide a weak dependence on the Kohn-Sham starting point and a much better agreement with reference calculations. The present findings rationalize the theory to experiment possible discrepancies observed in previous G{sub 0}W{sub 0} and Bethe-Salpeter studies of bulk water. The increase of the optical gap with increasing cluster size is consistent with the evolution from gas to dense ice or water phases and results from an enhanced screening of the electron-hole interaction.

37 INORGANIC, ORGANIC, PHYSICAL, AND ANALYTICAL CH↗

Electronic structure and metallization of cubic GdH{sub 3} under pressure: Ab initio many-body GW calculations

The electronic structures of the cubic GdH{sub 3} are extensively investigated using the ab initio many-body GW calculations treating the Gd 4f electrons either in the core (4f-core) or in the valence states (4f-val). Different degrees of quasiparticle (QP) self-consistent calculations with the different starting points are used to correct the failures of the GGA/GGA + U/HSE03 calculations. In the 4f-core case, GGA + G{sub 0}W{sub 0} calculations give a fundamental band gap of 1.72 eV, while GGA+ GW{sub 0} or GGA + GW calculations present a larger band gap. In the 4f-val case, the nonlocal exchange-correlation (xc) functional HSE03 can account much better for the strong localization of the 4f states than the semilocal or Hubbard U corrected xc functional in the Kohn–Sham equation. We show that the fundamental gap of the antiferromagnetic (AFM) or ferromagnetic (FM) GdH{sub 3} can be opened up by solving the QP equation with improved starting point of eigenvalues and wave functions given by HSE03. The HSE03 + G{sub 0}W{sub 0} calculations present a fundamental band gap of 2.73 eV in the AFM configuration, and the results of the corresponding GW{sub 0} and GW calculations are 2.89 and 3.03 eV, respectively. In general, for the cubic structure, the fundamental gap from G{sub 0}W{sub 0} calculations in the 4f-core case is the closest to the real result. By G{sub 0}W{sub 0} calculations in the 4f-core case, we find that H or Gd defects can strongly affect the band structure, especially the H defects. We explain the mechanism in terms of the possible electron correlation on the hydrogen site. Under compression, the insulator-to-metal transition in the cubic GdH{sub 3} occurs around 40 GPa, which might be a satisfied prediction.

71 CLASSICAL AND QUANTUM MECHANICS, GENERAL PHYSIC↗

Strongly hyperbolic quasilinear systems revisited, with applications to relativistic fluid dynamics

We revisit the theory of first-order quasilinear systems with diagonalizable principal part and only real eigenvalues, what is commonly referred to as strongly hyperbolic systems. We provide a self-contained and simple proof of local well-posedness, in the Hadamard sense, of the Cauchy problem. Our regularity assumptions are very minimal. As an application, we apply our results to systems of ideal and viscous relativistic fluids, where the theory of strongly hyperbolic equations has been systematically used to study several systems of physical interest.

Mathematics↗

Quantum computation of nanosheets in a background magnetic field for external control of nanosystems

We ran quantum simulations on lattices of interacting fermions in a background magnetic field [1]. We created discretized Hamiltonian matrices for each system and pre-processed them so that they could be read by the quantum simulator implemented in IBM's Qiskit package [2]. The simulator uses the variational quantum eigensolver (VQE) algorithm [3] to approximate the ground state energy of the input Hamiltonian, which is then compared to the classically-calculated exact value. This classical evaluation is possible because the systems studied are still within the grasp of classical computers. We also investigated ways to simplify the system before the pre-processing step. We found that in almost all cases these simplifications made our quantum calculations converge to better answers using simpler starting conditions. We are exploring other ways to interface with IBM's quantum computing software so that we no longer need to explicitly define the Hamiltonian matrices for our systems, a problem which will itself soon become intractable.

71 CLASSICAL AND QUANTUM MECHANICS, GENERAL PHYSIC↗

Prediction on Underestimation of Statistical Uncertainty in Monte Carlo Eigenvalue Calculation for Two-dimensional Heterogeneous Color Set Geometry

In a Monte Carlo calculation, statistical uncertainty is inevitably included in the result. The statistical uncertainty is commonly estimated by the central limit theorem without considering correlation among generations. The statistical uncertainty is usually expected to decrease at the inverse square root of the number of neutron histories. However, in a Monte Carlo eigenvalue calculation in a large geometry whose dominance ratio is close to unity, the central limit theorem would not be valid and thus the statistical uncertainty of fission source distribution can be significantly underestimated when inter-cycle correlation is not appropriately taken into account. The true statistical uncertainty can be calculated by statistical processing of the results obtained by multiple Monte Carlo runs using different initial random seeds, but this approach has a drawback of large calculation cost. Therefore, prediction of true statistical uncertainty by alternate approach is desirable. One of the solutions for this issue is an evaluation method of statistical uncertainty based on the eigenfunction expansion and the Autoregressive model (AR model). In this method, the true statistical uncertainty is theoretically derived by applying the AR model to behavior of expansion coefficients of eigenfunction in fission source distribution during cycles. Underestimation ratio of statistical uncertainty of fission source distribution can be predicted by this method. Note that the underestimation ratio of statistical uncertainty is defined by the ratio of estimated error by a Monte Carlo calculation to the true statistical error. In the previous study, the proposed method was verified in a simple condition such as one energy group, one dimensional homogeneous slab geometry. Contributions of this summary are as follows. Firstly, expression of the noise term that is appeared in the AR model is theoretically derived. Secondly, prediction on underestimation ratio of statistical uncertainty is carried out in a more general case, i.e., flux of each fuel cell in 2 x 2 heterogeneous assembly geometry taken from the C5G7 benchmark problem. Through the present study, applicability of the theoretical prediction for underestimation of statistical uncertainty in practical problems will be verified. (authors)

21 SPECIFIC NUCLEAR REACTORS AND ASSOCIATED PLANTS↗

On quark-antiquark approximation in light front QCD with zero gluon modes

We consider a transition to the light front Hamiltonian from theories quantized on spacelike planes approaching to the light front. In this approach we preserve the dynamics of zero mode present in the theories near the light front. We make the limit transition differently for zero and nonzero modes. This leads to the appearance of some phenomenological parameter which can be used to describe vacuum effects. Also we use a lattice gauge invariant regularization in transverse coordinate space. As an illustration of our scheme we consider the quark-antiquark bound states problem in 2+1 dimensions. We construct basis states in the light front Fock space and provide detailed computations of the Hamiltonian matrix in this basis. These steps allow us to construct the matrix eigenvalue equation. In conclusion we discuss the nuances of obtained results.

71 CLASSICAL AND QUANTUM MECHANICS, GENERAL PHYSIC↗

Sodium-Cooled Fast Reactor Reference Plant Model

This report details the progress of Idaho National Laboratory (INL) in creating a reference plant multiphysics model for the Advanced Burner Test Reactor (ABTR). This model was developed under Task 13 of the U.S. Nuclear Regulatory Commission project “Development and Modeling Support for Advanced Non-Light Water Reactors,” and is an extension of the reference plant model developed in Task 4b, which was improved upon in the following ways. (1) The discrete ordinates method was used in lieu of the super-homogenization (SPH)-corrected diffusion approximation in order to better capture the anisotropic scattering contribution and the neutron leakage change due to thermal expansion. (2) The novel neutronic spatial discretization approach, termed the ring-heterogeneous (RH) approximation, was conceptualized and introduced to capture the differential expansion of the materials in the core. This new technique proved capable of preserving fission rates and maintaining the eigenvalue within 2.5% and 266 pcm with 9 neutron energy groups, respectively. Separating the different materials in the core enables the differential expansion of materials to be explicitly accounted for, eliminating the need for problem-specific cross-section functionalization techniques. (3) The SAM model for the core and system thermal-hydraulics analysis was updated to include 61 channels instead of just four representative ones. This enables users to obtain improved spatial resolution for sodium temperature and density scalar fields. (4) All the mesh files were created via the Multiphysics Object-Oriented Simulation Environment (MOOSE) Reactor module, eliminating all reliance on external tools for mesh creation. (5) Finally, the fuel axial expansion now leverages the HT9 and UPuZr material properties that have been validated against experimental data. The reference plant model was used to perform a full-core unprotected loss of flow (ULOF) transient calculation, including neutronics, thermal and mechanical feedback mechanisms. Future work will be devoted to further enhancements of the model. Potential improvements to the model include the addition of the control rod driveline expansion feedback and the upgrading of the support plate model so as to explicitly include 3D effects. Additionally, a Nuclear Energy Advanced Modeling and Simulation funded parallel effort has completely automated the creation of the ring-heterogeneous (RH) mesh from the fully heterogeneous (FH) geometry, thus maximizing user friendliness for the sodium fast reactor sodium-cooled fast reactor (SFR) workflow and will be incorporated in future work.

22 GENERAL STUDIES OF NUCLEAR REACTORS↗

Data-driven modeling and control of dynamical systems using Koopman and Perron-Frobenius operators

This dissertation studies the data-driven modeling and control problem of nonlinear systems by exploiting the linear operator theoretic framework involving Koopman and Perro-Frobenius operator. A systematic linear-operator based controller design procedure has been established, which can be used to solve a variety of nonlinear control problems, including feedback stabilization using control Lyapunov functions, optimal quadratic regulation using Koopman eigenfunctions and convex optimization formulation of optimal control problem using P-F and Koopman operator approximation. As the core of data-driven modeling, we first propose a new algorithm for the finite-dimensional approximation of the linear transfer Koopman and Perron-Frobenius operator from time-series data. We argue that the existing approach for the finite-dimensional approximation of these transfer operators such as Dynamic Mode Decomposition (DMD) and Extended Dynamic Mode Decomposition (EDMD) do not capture two important properties of these operators, namely positivity and Markov property. The algorithm we propose preserves these two properties. We call the proposed algorithm as naturally structured DMD (NSDMD) since it retains the inherent properties of these operators. Naturally structured DMD algorithm leads to a better approximation of the steady-state dynamics of the system regarding computing Koopman and Perron- Frobenius operator eigenfunctions and eigenvalues. However, preserving positivity property is critical for capturing the real transient dynamics of the system. This positivity property of the transfer operators and it's finite-dimensional approximation play an important role for controller and estimator design of nonlinear systems. To solve the feedback stabilization problem for nonlinear control systems, we tried to take advantage of the Koopman operator framework. The Koopman operator approach provides a linear representation for a nonlinear dynamical system and a bilinear representation for a nonlinear control system. The problem of feedback stabilization of a nonlinear control system is then transformed to the stabilization of a bilinear control system. We propose a control Lyapunov function (CLF)-based approach for the design of stabilizing feedback controllers for the bilinear system. The search for finding a CLF for the bilinear control system is formulated as a convex optimization problem. This leads to a schematic procedure for designing CLF-based stabilizing feedback controllers for the bilinear system and hence the original nonlinear system. Another advantage of the proposed controller design approach outlined in this dissertation is that it does not require explicit knowledge of system dynamics. In particular, the bilinear representation of a nonlinear control system in the Koopman eigenfunction space can be obtained from time-series data. Next, we study the optimal quadratic regulation problem for nonlinear systems. The linear operator theoretic framework involving the Koopman operator is used to lift the dynamics of nonlinear control system to an infinite-dimensional bilinear system. The optimal quadratic regulation problem for nonlinear system is formulated in terms of the finite-dimensional approximation of the bilinear system. A convex optimization-based approach is proposed for solving the quadratic regulator problem for bilinear system. We applied a variety of examples and compared the simulation results between our framework and conventional LQR control using linearized model. For more general optimal control problems, we provide a density-function based convex formulation for the optimal control problem of the nonlinear system. The convex formulation relies on the duality result in the stability theory of a dynamical system involving density function and Perron-Frobenius operator. The optimal control problem is formulated as an infinite-dimensional convex optimization program. The finite-dimensional approximation of the optimization problem relies on the recent advances made in the data-driven computation of the Koopman operator, which is dual to the Perron-Frobenius operator. Simulation results are presented to demonstrate the application of the developed framework.

Huang, Bowen↗

A semi-analytical solver for the Grad-Shafranov equation

In toroidally confined plasmas, the Grad-Shafranov equation, in general, a non-linear partial differential equation, describes the hydromagnetic equilibrium of the system. This equation becomes linear when the kinetic pressure is proportional to the poloidal magnetic flux and the squared poloidal current is a quadratic function of it. In this work, the eigenvalue of the associated homogeneous equation is related with the safety factor on the magnetic axis, the plasma beta, and the Shafranov shift; then, the adjustable parameters of the particular solution are bounded through physical constrains. The poloidal magnetic flux becomes a linear superposition of independent solutions and its parameters are adjusted with a non-linear fitting algorithm. This method is used to find hydromagnetic equilibria with normal and reversed magnetic shear and defined values of the elongation, triangularity, aspect-ratio, and X-point(s). The resultant toroidal and poloidal beta, the safety factor at the 95% flux surface, and the plasma current are in agreement with usual experimental values for high beta discharges and the model can be used locally to describe reversed magnetic shear equilibria.

70 PLASMA PHYSICS AND FUSION TECHNOLOGY↗

Extended Applications of Subgrid Representation in the 2D/1D Method

Recent efforts in MPACT have focused on improving the performance of the 2D/1D subplane implementation to help target computational performance goals. Here, we build on previous efforts that targeted the use of subgrid treatments to improve the accuracy of control rod representation, presenting three additional applications of subgrid treatments with the goal of reducing the computational burden of simulations. These subgrid applications include treatment of spacer grids, thermal feedback, and axial reflector material representation. With these approaches, a single method of characteristics (MOC) plane can contain several different materials axially that are represented explicitly via subgrids on the coarse mesh finite difference (CMFD) mesh but are axially homogenized on the MOC mesh. This allows for a substantial reduction in the number of MOC planes needed in the calculation through the introduction of an approximate treatment, particularly with regard to the self-shielded cross sections and MOC-informed radial current coupling coefficients in CMFD. Several test problems ranging from single rod to quarter core are used to assess the solution accuracy and performance of these various subgrid representations. Overall, the accuracy of the approximations seems very reasonable, with extremely small differences in eigenvalue observed and maximum pin power errors in the 0.5% to 1.0% range. Several cases show substantial value in the compromise between accuracy and computational performance. Others highlight the new computational hurdles that future research will aim to resolve.

22 GENERAL STUDIES OF NUCLEAR REACTORS↗

Adapting CLUTCH methodology to multigroup TSUNAMI-3D for eigenvalue sensitivity calculations

The sensitivity of the eigenvalue to uncertainties in nuclear data and its evaluation are important for nuclear criticality safety. TSUNAMI-3D sequences within the SCALE code system offer several options to the user community for calculating eigenvalue sensitivity coefficients with multigroup (MG) and continuous energy (CE) 3D transport capabilities. TSUNAMI-3D sequences implement the adjoint-based perturbation theory with MG KENO code, the Contributon Linked eigenvalue sensitivity/Uncertainty estimation via Track length importance CHaracterization (CLUTCH) method with CE KENO code, and the Iterated Fission Probability (IFP) method with CE KENO and Shift codes. Each method has benefits and limitations depending on the problem that is run. The work presented here aims to adapt the CLUTCH method, which enables the Contributon method's mesh-free, memory-efficient approach for calculating adjoint-weighted tallies for sensitivity calculations, to the MG TSUNAMI-3D sequence. This application would eliminate the explicit adjoint KENO calculation, as well as the memory-consuming mesh flux moment tallies required by the conventional MG TSUNAMI-3D. Smaller memory footprints in the CLUTCH methodology and relatively shorter runtimes in MG KENO transport can make MG TSUNAMI-3D a viable method for some complex problems. Moreover, this adaptation allows MG sensitivity calculations with Shift, ORNL's next-generation high-performance Monte Carlo transport code, which currently does not offer any sensitivity capabilities with MG particle transport simulations. Initial implementation of the new MG TSUNAMI-3D sequence and its preliminary results with a selected critical benchmark experiment in the Verified, Archived Library of Inputs and Data (VALID) are presented in this study.

KENO↗

Rotational symmetry relation for efficient response function generation in the coarse mesh transport method COMET

The coarse mesh transport code COMET is a continuous energy hybrid stochastic-deterministic neutronics solver with high fidelity and formidable computation speed in solving reactor core problems. Its method is based on the incident flux expansion theory. In this work, we take advantage of the local geometric symmetry in many reactor cores lattices (e.g., fuel lattices and reflector blocks) to develop relations among the flux response expansion coefficients for symmetric surfaces to further improve the computational efficiency of the COMET response function generation tool (method). This is done by a rigorous derivation of the transformation matrices for the angular and spatial expansion moments resulting from a rotation of a coarse mesh by an arbitrary angle. The relations for the response coefficients for the symmetric surfaces can be then written as the Kronecker product of those transformation matrices. The method is implemented into COMET and tested on two advanced high temperature reactor (AHTR) full-length single assembly benchmark problems. The COMET results using the response function library based on the symmetry relations were compared to those using the library directly generated by continuous energy Monte Carlo for all surfaces. It was found that the eigenvalues and stripe-wise fission densities using the two libraries are in statistical agreement as expected. This indicates that the new method maintains the high fidelity of the original COMET method while improving the computational efficiency in the response function generation by 270% to 400%, depending on the local geometric symmetry. This method also reduces the size of the response function library by the same magnitude (270% to 400%). (authors)

21 SPECIFIC NUCLEAR REACTORS AND ASSOCIATED PLANTS↗

Year End Progress Report on Rattlesnake Improvements

Rattlesnake is a MOOSE-based radiation transport application developed at INL to support modern multi-physics simulations. At the beginning of the last year, Rattlesnake was able to perform steady-state, transient and eigenvalue calculations for the multigroup radiation transport equations. Various discretization schemes, including continuous finite element method (FEM) with discrete ordinates method (SN) and spherical harmonics expansion method (PN) for the self-adjoint angular flux (SAAF) formulation, continuous FEM (CFEM) with SN for the least square (LS) formulation, diffusion approximation with CFEM and discontinuous FEM (DFEM), have been implemented. A separate toolkit, YAKXS, for multigroup cross section management was developed to support Rattlesnake calculations with feedback both from changes in the field variables, such as fuel temperature, coolant density, and etc., and in isotope inventory. The framework for doing nonlinear diffusion acceleration (NDA) within Rattlesnake has been set up, and both NDA calculations with SAAF-SN-CFEM scheme and Monte Carlo with OpenMC have been performed. It was also used for coupling BISON and RELAP-7 for the full-core multiphysics simulations. Within the last fiscal year, significant improvements have been made in Rattlesnake. Rattlesnake development was migrated into our internal GITLAB development environment at the end of year 2014. Since then total 369 merge requests has been accepted into Rattlesnake. It is noted that the MOOSE framework that Rattlesnake is based on is under continuous developments. Improvements made in MOOSE can improve the Rattlesnake. It is acknowledged that MOOSE developers spent efforts on patching Rattlesnake for the improvements made on the framework side. This report will not cover the code restructuring for better readability and modularity and documentation improvements, which we have spent tremendous effort on. It only details some of improvements in the following sections.

22 GENERAL STUDIES OF NUCLEAR REACTORS↗

Development of a Library for Conducting Monte Carlo Tallies on Heterogeneous Systems

A library is introduced that enables users to conduct Monte Carlo tallies on heterogeneous computing systems. The library can be built with or linked to a Monte Carlo transport code to provide tally capabilities for the project. The library provides a set of tools and an interface for constructing tallies that can be computed on CPUs or GPUs regardless of where the sample data was generated. The library can be incorporated into existing Monte Carlo transport codes at both the research and production level with minimal intrusion into the transport code. The tally library is capable of computing a variety of tallies including track-length histogram (mesh) tallies, collisionbased kernel density estimators, functional expansion tallies, eigenvalue sensitivity coefficients via iterated fission probability, and sensitivities of reaction rates or adjoint-weighted quantities via differential operator sampling. The parallelization methods are discussed and speedups are shown for 2-D and 3-D problems with comparisons to both the tally library CPU implementations as well as the tally capabilities in the host transport code. Performance results and examples are shown for MCNP, MCATK, and OpenMC.

97 MATHEMATICS AND COMPUTING↗

Analyzing Prospects for Quantum Advantage in Topological Data Analysis

Lloyd [Nat. Commun. , 10138 (2016)] were first to demonstrate the promise of quantum algorithms for computing Betti numbers, a way to characterize topological features of data sets. Here, we propose, analyze, and optimize an improved quantum algorithm for topological data analysis (TDA) with reduced scaling, including a method for preparing Dicke states based on inequality testing, a more efficient amplitude estimation algorithm using Kaiser windows, and an optimal implementation of eigenvalue projectors based on Chebyshev polynomials. We compile our approach to a fault-tolerant gate set and estimate constant factors in the Toffoli complexity. Our analysis reveals that superquadratic quantum speedups are only possible for this problem when targeting a multiplicative error approximation and the Betti number grows asymptotically. Further, we propose a dequantization of the quantum TDA algorithm that shows that having exponentially large dimension and Betti number are necessary, but insufficient conditions, for superpolynomial advantage. We then introduce and analyze specific problem examples which have parameters in the regime where superpolynomial advantages may be achieved, and argue that quantum circuits with tens of billions of Toffoli gates can solve seemingly classically intractable instances. Published by the American Physical Society 2024

97 MATHEMATICS AND COMPUTING↗

Better band gaps with asymptotically corrected local exchange potentials

In this study, we formulate a spin-polarized van Leeuwen and Baerends (vLB) correction to the local density approximation (LDA) exchange potential [R. van Leeuwen and E. J. Baerends, Phys. Rev. A 49, 2421 (1994)] that enforces the ionization potential (IP) theorem following T. Stein et al. [Phys. Rev. Lett. 105, 266802 (2010)]. For electronic-structure problems, the vLB correction replicates the behavior of exact-exchange potentials, with improved scaling and well-behaved asymptotics, but with the computational cost of semilocal functionals. The vLB + IP correction produces a large improvement in the eigenvalues over those from the LDA due to correct asymptotic behavior and atomic shell structures, as shown in rare-gas, alkaline-earth, zinc-based oxides, alkali halides, sulfides, and nitrides. In half-Heusler alloys, this asymptotically corrected LDA reproduces the spin-polarized properties correctly, including magnetism and half-metallicity. We also consider finite-sized systems [e.g., ringed boron nitride (B 12 N 12 ) and graphene (C 24 )] to emphasize the wide applicability of the method.

72 PHYSICS OF ELEMENTARY PARTICLES AND FIELDS↗

A mixed, unified forward/inverse framework for earthquake problems: fault implementation and coseismic slip estimate

SUMMARY We introduce a new finite-element (FE) based computational framework to solve forward and inverse elastic deformation problems for earthquake faulting via the adjoint method. Based on two advanced computational libraries, FEniCS and hIPPYlib for the forward and inverse problems, respectively, this framework is flexible, transparent and easily extensible. We represent a fault discontinuity through a mixed FE elasticity formulation, which approximates the stress with higher order accuracy and exposes the prescribed slip explicitly in the variational form without using conventional split node and decomposition discrete approaches. This also allows the first order optimality condition, that is the vanishing of the gradient, to be expressed in continuous form, which leads to consistent discretizations of all field variables, including the slip. We show comparisons with the standard, pure displacement formulation and a model containing an in-plane mode II crack, whose slip is prescribed via the split node technique. We demonstrate the potential of this new computational framework by performing a linear coseismic slip inversion through adjoint-based optimization methods, without requiring computation of elastic Green’s functions. Specifically, we consider a penalized least squares formulation, which in a Bayesian setting—under the assumption of Gaussian noise and prior—reflects the negative log of the posterior distribution. The comparison of the inversion results with a standard, linear inverse theory approach based on Okada’s solutions shows analogous results. Preliminary uncertainties are estimated via eigenvalue analysis of the Hessian of the penalized least squares objective function. Our implementation is fully open-source and Jupyter notebooks to reproduce our results are provided. The extension to a fully Bayesian framework for detailed uncertainty quantification and non-linear inversions, including for heterogeneous media earthquake problems, will be analysed in a forthcoming paper.

58 GEOSCIENCES↗