Casper's Research Plan 2023

# Casper's Research Plan 2023/24 > This document lists my plans for late 2023/early 2024, a timeline and a to-do list, as well my responsibilities, and perhaps some meeting notes. Anyone with the link can edit! ## Table of Contents [ToC] ## Projects/Papers ### StaggeredKernels.jl: data structures and utility functions for operating on fields on staggered grids #### Repository https://github.com/cpranger/StaggeredKernels.jl #### Status Functional, but not yet portable/performant #### To do list - [x] Publish StaggeredKernels.jl to a public repository [here](https://github.com/cpranger/StaggeredKernels.jl) - [ ] Implement [tinykernels.jl](https://github.com/utkinis/TinyKernels.jl) as a backend to StaggeredKernels.jl, and test in the real-world on GPU ### DamageBreakage.jl: algorithms for solving time-dependent, coupled, nonliear parabolic/hyperbolic equations matrix-free with robust error control and convergence behavior #### Repository - https://github.com/cpranger/DamageBreakage #### Status Getting there, implementing time integration schemes now. #### To do list - [x] Publish Damage Breakage research to public repository [here](https://github.com/cpranger/DamageBreakage) - [x] Test the vector identities that ought to be conserved on the staggered grid [here](https://github.com/cpranger/DamageBreakage/blob/main/scripts/test-de-rham.jl) - [x] Test the analytically computed eigenmodes for scalar Laplacian and the two vector Laplacians for curl-free and divergence-free flows [here](https://github.com/cpranger/DamageBreakage/blob/main/scripts/test-eigenmodes.jl) - [x] Implement power iteration and assert that the power iteration algorithm yields the same results numerically [here](https://github.com/cpranger/DamageBreakage/blob/main/scripts/test-powerit.jl) - [x] Implement and test a Chebyshev solver, which relies on the spectral bounds given by the power iteration algorithm [here](https://github.com/cpranger/DamageBreakage/blob/main/scripts/test-chebyshev.jl). Conclusion: power iteration takes at least as much computational work as the Chebyshev iterations themselves, and introduces a hell of a lot of inner products, which some people are allergic to - [x] Implement and test Lanczos iteration [here](https://github.com/cpranger/DamageBreakage/blob/main/scripts/test-lanczos.jl). Conclusion: extremely elegant method, but progressively approaches the spectral bounds from the inside, which makes Chebyshev unstable - [x] Implement and test Rayleigh quotient iteration. Conclusion: needs an iterative solver like Chebyshev to work in the first place, and recursing like that just leads to instability - [x] Merge and test Conjugate Gradients (CG) with Lanczos iteration [here](https://github.com/cpranger/DamageBreakage/blob/main/scripts/test-cg.jl). Conclusion: at the cost of a few inner products, one now doesn't need to know the operator's spectral bounds, but will get it nevertheless for free for use elsewhere (e.g. error estimation, time stepping) - [x] Implement and test linearization by complex step differentiation [here](https://github.com/cpranger/DamageBreakage/blob/main/scripts/test-linearization.jl) - [x] Implement and test Schur complement reduction here [here](https://github.com/cpranger/DamageBreakage/blob/main/scripts/test-churcomplement.jl) - [x] Implement and test the TR-BDF2 time integration scheme with parabolic and hyperbolic scalar and vector problems [here](https://github.com/cpranger/DamageBreakage/blob/main/scripts/test-integration.jl) - [x] Check PETSc KSP source for CG -- eigenvalues are optionally computed there - [x] Check out symmetric preconditioning by the diagonal - [ ] Implement and test simple absorbing BC. - [ ] Couple a hyperbolic and a parabolic problem using the developed IMEX scheme - [ ] Construct 1D/2D toy model - [x] Fix implementation of constitutive tensor - [ ] ~~Obtain the previous by complex step differentiation? (e.g. $\rho \frac{dv}{dt} = \nabla \cdot \left(\frac{\partial U}{\partial e} + \frac{\partial U}{\partial \nabla \alpha} \right)$)~~ ### Technical paper for geodynamics people on te latter #### About I've formed a lot of ideas on how to solve the kind of systems of equations that are of interest to geodynamicists in an iterative, matrix-free way that is popular right now. It would be good to write those down. #### Repository Will be provided. submit to G-cubed #### Figures/examples - [ ] Intermediate-depth earthquakes with Arne/Marcel - [ ] Damage-Breakage #### Status Quite a lot has been written down, but is spread out over multiple paper drats, pdfs, mathematica note books, hand-written notes. #### To do list - [ ] Make to do list ### Final report TEAR #### About #### Repository #### Status #### To do list - [ ] Make to do list ### Collaboration in with Marcel and Arne on intermediate-depth earthquakes #### Repository #### Status Arne and Marcel agree to collaborate on intermediate-depth earthquakes, making it possible that progress on Damage-Breakage translates into progress on intermediate-depth earthquakes before the end of the TEAR project. #### To do list - [x] Contact them, ask how it's going - [ ] Check in on Marcel's progress with Newton solver + PT some time after GEOMOD - [ ] Provide Arne and Marcel with a simple 1D toy model of diffusion-regularized thermal runaway coupled to the elastic wave equation - [ ] (Boundary conditions for 1D model? --> corner case that can lead to indefinite friction drop) ## Last 2 weeks' work - CG algorithm compared in detail against PETSc implementation, are largely equivalent and produce equivalent results. Translated PETSc implementation not comitted to repository, but saved elsewhere. - Systems made solvable by CG by considering the following: - CG converges fine under appropriately scaled Neumann and Dirichlet conditions, provided the initial guess satisfies the latter - Out of the states (velocity, elastic strain, rheological params α & β), Neumann conditions are appropriate for the latter two, boundary conditions are not needed on elastic strain (due to staggering and because eliminated by Schur), and both Dirichlet and Neumann conditions can be appropriate for the velocity field - Boundary conditions can be applied in rate form directly to $F$ in $\dot{y} = F(y)$, rather than implicitly on the state $y$ by solving the minimization problem |R = A(x) - B| = 0. This makes initial guessing the correct boundary solution trivial - variable loading and simple absorbing boundary conditions can still conveniently be applied through the RHS on velocity - IMEX TR-BDF2 time integration with embedded error computation running correctly except on hyperbolic problems, which are solved by Schur complement reduction and need just a bit more boilerplate code (i.e. to be able to compute TR-BDF2 stages on objects of type SchurComplement -- to be finalized tomorrow). - Schur complement reduction is independently verified to function correctly.