Fusion

Fusion Personnel: Steve Jardin (PPPL), Hank Strauss (PPPL), Jin Chen (PPPL)

The goal of this effort is to develop high-order adaptive discretization technologies aimed toward solving MHD problems of importance to the fusion community. Efforts to date have involved examining the M3D code used at PPPL with an eye toward its installation at SCOREC and implementing a specific MHD capability within the SCOREC (RPI) Trellis framework to provide an initial demonstration of the type of procedures that can be developed.

The first activity included examining the M3D code for the purposes of (i) gaining a better understanding of the class of MHD problem being solved by PPPL, (ii) gaining an initial understanding of the numerical methods used in M3D and the methods of implementation used, (iii) determining the level of effort required to install M3D on the SCOREC computer system. An important lesson learned during this process is that the class of problems solved by M3D requires the application of advanced methods and that because of the evolutionary implementation of these methods, there is a complex interaction of the numerical methods, discretization technologies, and (parallel) computing environment. To effectively implement high-order and adaptive methods into M3D requires a stronger separation of the discretization library operations and the specific numerical methods used to solve the discretization. Due to the complexity of the class of problem being solved, it is not known in advance how easily or effectively those tasks can be undertaken.

As the first concrete step toward our long term goal, the SCOREC team has implemented one specific set of MHD equations previously considered by Hank Strauss and D. Longcope (J. Comp. Physics, 1998, 147) in the Trellis simulation system. Trellis is a geometry-based system for the implementation of high-order adaptive methods which already includes various high-order basis functions and allows the effective implementation of new weak forms and their discretization. The code is modular and individual components developed as part of this effort be transferred to and re-used in other simulation codes.

Using Trellis, we have investigated the two-dimensional, incompressible MHD equations in both the primitive variable formulation and stream function formulation. In our incompressible MHD model, both magnetic and velocity fields are solenoidal (divergence free) and there are two ways to enforce this constraint; (i) from the interior, using potentials, and (ii) from the exterior, using Lagrange multipliers.

We first used Lagrange multipliers to impose the incompressibility constraints because the magnetic part of the MHD model was trivial to add to our existing incompressible Navier-Stokes implementation in Trellis. We used semi-implicit time discretization and compute convective terms explicitly so that the stability in time is constrained by a CFL condition based on the fluid velocity. We used 2nd and 3rd order finite elements to discretize the primary fields (velocity and magnetic field) and spaces that are compatible with the Babuska-Brezzi conditions for the Lagrange multipliers. We were able to produce stable results on different grids and this first experiment enabled us to familiarize ourselves with MHD flows phenomenology.

In the stream function formulation, the equations consist of two Poisson equations and two time dependent nonlinear equations, one purely advective and the other advective diffusive. The nonlinearity is present in the form of Poisson brackets between two variables. In the weak form, the two time dependent equations are stabilized using the Streamline Upwind Petrov Galerkin (SUPG) method. Finite elements based on the Lagrange shape functions are used to discretize the weak form in space and a Backward Euler scheme is used for the time integration. The solution procedure consists of solving one equation at a time. This allows us to linearize the nonlinear terms in the equations by using values at the previous time steps for one of the variables in the Poisson bracket.

Our test case is the tilt instability example (see figure below) studied by Strauss and Longcope for the incompressible case and by Richard, Sydora and Ashour-Abdalla for the compressible case. The initial equilibrium state corresponds to a bipolar vortex for the magnetic field with zero initial velocity for the fluid. Two oppositely directed currents are embedded in a constant magnetic field in the horizontal plane. As a result, there are two sets of closed magnetic field line loops or islands pressed to each other in equilibrium by the action of the magnetic field. This equilibrium is however unstable under small perturbations. In this case, the two islands rotate until reaching a horizontal position from which they are expelled in opposite directions under the magnetic field influence.