“Assimilating” Solar Data into MHD Models of the Solar Atmosphere

1. “Assimilating” Solar Data into MHD Models of the Solar Atmosphere W.P. Abbett SSL UC Berkeley HMI Team Meeting, Jan 2005

2. A Quick Synopsis of what’s needed: An MHD model active region requires both • A means of specifying the photospheric boundary in a way consistent with the observed evolution of the vector magnetic field at the photosphere. • An initial atmosphere throughout the computational volume, which must encompass the β~1 plasma of the photosphere along with the low-β coronal plasma. • A means of specifying the other external boundaries.

3. What we Have: • With new means of specifying physically consistent photospheric flows (ILCT, MEF, MSR), we have u and B (and thus E┴) at z=0 for all t. • Using the Wheatland technique, we can generate a divergence-free, “force-free” magnetic field throughout the volume that matches the imposed boundary conditions at t = 0.

4. What we still need: • Inversion techniques like ILCT and MEF do not provide the vertical gradients of u (since a single vector magnetogram gives no information about the vertical gradients of B) nor do they say anything about the evolution of the transverse components of the magnetic field. • Yet higher than first-order accurate numerical schemes require information from cells below the boundary layer to properly calculate the fluxes at the boundary face necessary to update the solution.

5. The Strategy:

6. The “Active” Boundary • Treat as a code-coupling exercise --- we now have two distinct 3D regions: • An MHD model with a domain that encompasses the layers above the photosphere and extends out to the low corona, and • A 3D dynamic “boundary” layer, with its lower boundary at the photosphere, and its upper boundary at the base of the MHD model corona

7. The “Active” Boundary • Make a physical assumption: Coronal forces do not affect photospheric motions. • Then in the “boundary” layer, we can assume that ILCT or MEF flows permeate the entire layer, and implicitly solve (using the ADI method) the induction continuity, and (simplified energy) equations given the prescribed flow field.

8. Coupling the codes • Use a modified version of this ADI “boundary” code to implicitly solve the induction, continuity, and (a simple) energy equation given the ILCT or MEF prescribed flow field • The upper boundary of the ADI code extends into the lower active zones of the MHD model, and the ghost cells of the MHD model are specified by the upper active zones of the ADI code. • The cadence is determined by the (explicit) MHD code.

9. Timestep severely restricted by CFL condition up in the corona --- the characteristic flow speed in the corona far exceeds that of the photosphere. • Simple scaling of e.g. the update cadence, field strength, or size scale of AR will not suffice, if one desires to maintain the physics of the transition layers ∂B/∂t =  x (v x B)

10. Principal Challenges: • Extreme separation of time scales CFL condition: Δt < Δz/c Sound speeds large in the low corona; Characteristic flow speeds along loops can be 100’s of kms per second • Extreme separation of spatial scales Evolution of the AR’s magnetic field evolution is not independent of the global magnetic field (8210 connected by a trans-equatorial loop to another AR complex!) Convective granulation pattern small compared to the size scale of AR at the photosphere

11. Methods of Attack: • Explicit temporal differencing • Solve for qn+1 directly from the state of the atmosphere at qn • Fully implicit differencing • Define qn+1 implicitly in terms of qn --- requires the inversion of a large, sparse matrix • Semi-implicit techniques

12. Advantages/Disadvantages of each approach • Explicit methods: Accurate, but numerically stable only if CFL condition is satisfied • Semi-implicit methods: efficient when “stiffness” comes from linearities. • Coronal plasma • Fully-implicit methods: efficient when “stiffness” comes from non-linearities • Timestep restricted by the dynamic timescale of interest

13. Fully-implicit technique: NR • Widely used in complex non-linear 1D problems (e.g., non-LTE radiation hydro) • Basic idea: Multi-dimensional Taylor expansion of F(q) about the current state vector qk= (ρ, px, py, pz, ε, Bx, By, Bz): F(qk+1) = F(qk) + (dF/dq)|k (qk+1 - qk) + … • Then for 2nd order corrections solve: (dF/dq)|kδqk = - F(qk) then update qk : qk+1 = qk + δqk

14. Why is this not in general use for 3D MHD problems? • The calculation and storage of J = (dF/dq)|k is computationally prohibitive! • J is an N X N matrix where N = neq∙nx∙ny∙nz But there is a solution to this problem!! • “Jacobian-free” Newton-Krylov methods (new technique used in fluid dynamics, Fokker-Planck codes, and even 2D Hall MHD)

15. The Basic Idea: • Make an initial guess for the correction vector δq and form a “linear residual” r0: • Solve for δq using a “Krylov-based GMRES” technique: where βi is given by the minimization of in a least squares sense. The Important Point: We need only to calculate a matrix vector product Jv which we approximate using

16. Caveats • Can’t get something for nothing! The GMRES convergence rate slows as the timestep is increased • Solution: “pre-condition” δq before Krylov iteration. • More info about technique: see Knoll & Keyes JcP 2004 193,57.

17. Initial Tests: • IMHD (Implicit MHD) • F95 3D resistive finite-volume MHD code • Not “operator split” • Set up to interface with PARAMESH (to address the problem of disparate spatial scales) • Efficient use of memory (can run on a laptop) • Can generate static solutions to improve initial configuration • In it’s infancy! “First light” only a few weeks ago…..

18. Some tests…..