Jacobian-free Newton-Krylov Method Research Articles

Jacobian-free Newton-Krylov method for multilevel nonlocal thermal equilibrium radiative transfer problems

Context. The calculation of the emerging radiation from a model atmosphere requires knowledge of the emissivity and absorption coefficients, which are proportional to the atomic level population densities of the levels involved in each transition. Due to the intricate interdependence of the radiation field and the physical state of the atoms, iterative methods are required in order to calculate the atomic level population densities. A variety of different methods have been proposed to solve this problem, which is known as the nonlocal thermodynamical equilibrium (NLTE) problem. Aims. Our goal is to develop an efficient and rapidly converging method to solve the NLTE problem under the assumption of statistical equilibrium. In particular, we explore whether the Jacobian-Free Newton-Krylov (JFNK) method can be used. This method does not require an explicit construction of the Jacobian matrix because it estimates the new correction with the Krylov-subspace method. Methods. We implemented an NLTE radiative transfer code with overlapping bound-bound and bound-free transitions. This solved the statistical equilibrium equations using a JFNK method, assuming a depth-stratified plane-parallel atmosphere. As a reference, we also implemented the Rybicki & Hummer (1992) method based on linearization and operator splitting. Results. Our tests with the Fontenla, Avrett and Loeser C model atmosphere (FAL-C) and two different six-level Ca II and H I atoms show that the JFNK method can converge faster than our reference case by up to a factor 2. This number is evaluated in terms of the total number of evaluations of the formal solution of the radiative transfer equation for all frequencies and directions. This method can also reach a lower residual error compared to the reference case. Conclusions. The JFNK method we developed poses a new alternative to solving the NLTE problem. Because it is not based on operator splitting with a local approximate operator, it can improve the convergence of the NLTE problem in highly scattering cases. One major advantage of this method is that it is expected to allow for a direct implementation of more complex problems, such as overlapping transitions from different active atoms, charge conservation, or a more efficient treatment of partial redistribution, without having to explicitly linearize the equations.

An unconditionally stable scheme for the immersed boundary method with application in cardiac mechanics

The stability of the immersed boundary (IB) method is a challenge in simulating fluid–structure interaction problems, where time step constraints are significantly stricter than in pure fluid simulations. We propose a novel unconditionally stable scheme for the immersed boundary finite element (IBFE) method. The structure is handled implicitly and characterized by strain energy functions, rather than being modeled as fibers or membranes. Through energy estimate, we prove the unconditional stability of the fully discrete approximation in the absence of the convective term. In real simulations of cardiac mechanics problems, the time step is much larger, only limited by the Courant–Friedrichs–Lewy condition of the fluid. The novelty of this work lies in the combination of dual interpolation and distribution operators, the Jacobian-free Newton–Krylov method for solving nonlinear algebraic systems, and the semi-Lagrangian method for handling the convective term. To validate the effectiveness and accuracy of our approach, we present various benchmarks and conduct a quasi-static simulation of a three-dimensional real left ventricular model. We have shown that the numerical stability of our scheme is very robust even with much larger time step compared to conventional explicit IB methods. Our work paves the way for further works on efficient solvers of the IBFE method.

A Parallel Coloring Newton-Krylov Method for Multiphysics Coupling System in Nuclear Reactors

The Newton-Krylov method with the explicit Jacobian matrix is an efficient numerical method for solving the nuclear reactor nonlinear multiphysics coupling system. Compared with the Jacobian-free Newton-Krylov (JFNK) method, it has a better preconditioner matrix (the Jacobian matrix itself) and can achieve a more stable and faster convergence. How to compute the Jacobian matrix efficiently is a key issue for this method. The graph coloring algorithm is an essential technique and has been used to reduce the Jacobian computational burden by exploiting its sparsity. The fewer the coloring numbers in the Jacobian, the less the Jacobian computational cost will be. Besides, when computing the Jacobian in a distributed memory parallel environment, the parallel graph coloring algorithms are required because the Jacobian is distributed among processors. Currently, a popular parallel graph coloring algorithm has been used to color the Jacobian. However, this parallel graph coloring algorithm shows poor scalability in parallel. The coloring numbers will increase with the processors, resulting in poor Jacobian computational efficiency. In this paper, a more efficient parallel graph coloring method is proposed that aims to reduce the coloring numbers and improve Jacobian computation efficiency in parallel. The main feature of the new method is that the coloring numbers decrease with the increasing number of processors. A neutronics/thermal-hydraulic coupling problem arising from the simplified high-temperature gas coolant model is utilized to assess the performance of the newly proposed method. The results show that (1) the parallel coloring number is reduced significantly, (2) the Jacobian computed by the new method is completely correct and excellent parallel scalability is achieved, and (3) the parallel coloring Newton-Krylov method with explicit Jacobian is more efficient and more stable than the parallel JFNK due to a better preconditioner.

Multi-Passage harmonic balance method for flutter prediction in multi-row configurations

A multi-passage harmonic balance method is presented for flutter prediction in multi-row configurations, and particular focus is on simulation of more than two blade rows. The unsteady flows in such configurations contain two kinds of multiple perturbations: one refers to those who have different frequencies and the other one refers to those who have the same frequency but different inter-blade phase angles (IBPAs). Most of current harmonic balance methods only consider the former one and the applications have been restricted to two-blade-row configurations. In the few studies that consider more rows, a one-to-one correspondence assumption between frequency and IBPA is usually made putting extra limitation on the used modes. This motivates our work, and a simple but effective idea is proposed by adding passages adopted per row. The theoretical basis of this work is provided by the spinning mode analysis. The present method is developed based on the time-spectral and multi-frequential form of harmonic balance. The treatment of pitchwise and inter-row boundary conditions in the case of multi passages are particularly addressed. A recently-developed Jacobian-free Newton-Krylov method with a simplified LU-SGS preconditioner is employed to solve the resultant harmonic balance equations. For validation, two typical multi-row flutter cases are investigated. The numerical results show the validity and accuracy of the developed method.

On the Use of the Jacobian-Free Newton Krylov Method to Generate One-Group Discontinuity and Super Homogenization Factors for Full-Core Neutron Diffusion Simulations

Coupled Monte Carlo (MC) and thermal-hydraulic analysis is valuable as a design or reference tool but can be slow, especially when implemented in a Picard iteration. Previous work has developed a novel prediction block to achieve convergence with fewer MC simulations. The prediction works in two stages: (1) a surrogate-like model predicts macroscopic cross sections on the fly and (2) a reduced-order neutronic model predicts the flux response to the updated cross sections. The main challenge with the prediction block is that the reduced-order neutronic model cannot reproduce the spatial flux distribution with high fidelity. This paper investigates the well-established Jacobian-free Newton Krylov (JFNK) method to preserve equivalency between a homogeneous (nodal diffusion) solution and a high-fidelity transport (MC) solution. Instead of performing multiple computationally consuming MC simulations, the nonlinear iterative approach iterates on correction parameters, e.g., assembly discontinuity factors (ADFs) or super homogenization (SPH) factors, using unexpensive nodal solutions. The JFNK approach does not require additional overhead from the MC solver to generate flux tallies. Further, the approach iterates on diffusion solutions produced directly from a desired code, thus ensuring that the parameters are compatible with that code. The approach is applied to correct a nodal diffusion solution for a realistic three-dimensional pressurized water reactor core. The results obtained in the paper show that the method is very successful in reproducing the heterogeneous solution (up to 2.5% difference in assembly flux for SPH and 0.3% for ADFs) without needing to modify the source code of the nodal diffusion solver. In addition, the results show that ADFs yield the best agreement and are also stable (i.e., weakly varying) when thermal-hydraulic fields are perturbed.

An Innovational Jacobian-Split Newton-Krylov Method Combining the Advantages of the Jacobian-Free Newton-Krylov Method and the Finite Difference Jacobian-Based Newton-Krylov Method

The Jacobian-free Newton-Krylov (JFNK) method is a widely used and flexible numerical method for solving the neutronic/thermal-hydraulic coupling system. The main property of JFNK is that the Jacobian-vector product is evaluated approximately by finite difference, avoiding the forming and storage of Jacobian explicitly. However, the lack of an efficient preconditioner is a major bottleneck for the JFNK method, leading to poor convergence. The finite difference Jacobian-based Newton-Krylov (DJNK) method is another advanced numerical method, in which the Jacobian matrix is formed and stored explicitly. The DJNK method can provide a better preconditioner for Krylov iteration than JFNK. However, how to compute the Jacobian matrix efficiently and automatically is a key issue for the DJNK method. By fully utilizing the sparsity of the Jacobian matrix and graph coloring algorithm, the Jacobian can be computed efficiently. Unfortunately, when there are dense rows/blocks, a huge computational burden will emerge due to the lack of sparsity, resulting in the extremely poor efficiency of Jacobian computation. In this work, a Jacobian-split Newton-Krylov (JSNK) method is proposed to resolve the dense row/block problem by combining the advantages of JFNK and DJNK. The main feature of the JSNK method is to split the Jacobian matrix into sparse and dense parts. The sparse part of the Jacobian matrix is explicitly constructed using the graph coloring algorithm while for the dense part, the Jacobian-vector product is approximated by finite difference. The computational complexity of the JSNK method is analyzed and compared to the JFNK method and the DJNK method from theoretical and experimental aspects and under different meshes. A simplified two-dimensional (2-D) high-temperature gas-cooled reactor (HTR) model and a simplified 2-D pressurized water reactor model are utilized to demonstrate the superiority of the JSNK method. The numerical results show that the JSNK method successfully resolved the dense rows/blocks. More importantly, its efficiency significantly outperforms the JFNK method and the DJNK method.

A GPU-Accelerated Method for 3D Nonlinear Kelvin Ship Wake Patterns Simulation

The study of ship waves is important for ship detection, coastal erosion and wave drag. This paper proposed a highly paralleled numerical computation method for efficiently simulating three-dimensional nonlinear kelvin waves. First, a numerical model for nonlinear ship waves is established based on potential flow theory, the Jacobian-free Newton–Krylov (JFNK) method and the boundary integral method. To reduce the amount of data stored in the JFNK method and improve the computational efficiency, a banded preconditioner method is then developed by formulating the optimal bandwidth selection rule. After that, a Graphics Process Unit (GPU)-based parallel computing framework is designed, and we used the Compute Unified Device Architecture (CUDA) language to develop a GPU solution. Finally, numerical simulations of 3D nonlinear ship waves under multiple scales are performed by using the GPU and CPU solvers. Simulation results show that the proposed GPU solver is more efficient than the CPU solver with the same accuracy. More than 66% GPU memory can be saved, and the computational speed can be accelerated up to 20 times. Hence, the computation time for Kelvin ship waves simulation can be significantly reduced by applying the GPU parallel numerical scheme, which lays a solid foundation for practical ocean engineering.

Development of a Fully-Implicit Second-Order system Thermal-Hydraulics analysis code for reactor modeling and simulations

There has been continuous advancement in numerical techniques and software engineering in the field of system thermal-hydraulics over the past decade. The current state-of-the-art involves the use of fully-implicit high-resolution schemes in conjunction with the Jacobian-Free Newton-Krylov (JFNK) method for both single-phase and two-phase flows. However, the application of these advanced schemes in system thermal-hydraulics codes is not a common practice. Given the opportunity and the need to develop a new modern system thermal–hydraulic code from scratch, there is no doubt that these recent advancements should be leveraged. In this work, we introduce a new modern object-oriented system thermal-hydraulics code named RETA. RETA is characterized by several key features, including: second-order temporal and spatial discretization schemes for both incompressible and compressible flows, as well as single-phase and two-phase flows; fully-implicit solution schemes utilizing both the Newton and PJFNK methods; and an object-oriented design of major fluid components and physical models, which allows for easy extension. A series of verification and demonstration studies are conducted to verify the implementation, demonstrate the accuracy improvement, and evaluate the performance of the new code. It is concluded that results from the second-order scheme are generally more reliable. Without the concern of a code stability issue, the second-order scheme should be the preferred and default choice for a code user or an analyst.

An implicit particle code with exact energy and charge conservation for electromagnetic studies of dense plasmas

A collisional particle code based on implicit energy- and charge-conserving methods is presented. A modified version of the particle-suppressed Jacobian-Free Newton-Krylov method that can enhance the solver efficiency is introduced. Mathematically, it is shown that this new approach can be viewed as a fixed-point iteration method for the particle positions. The model can exactly conserve global energy and local charge and can efficiently use time steps larger than the plasma period. The algorithm's ability to simulate dense plasmas accurately and efficiently is quantified by simulating the dynamic compression of a plasma slab via a magnetic piston in 1D planar geometry.

A Jacobian-Free Newton–Krylov Method to Solve Tumor Growth Problems with Effective Preconditioning Strategies

A Jacobian-free Newton–Krylov (JFNK) method with effective preconditioning strategies is introduced to solve a diffusion-based tumor growth model, also known as the Fisher–Kolmogorov partial differential equation (PDE). The time discretization of the PDE is based on the backward Euler and the Crank–Nicolson methods. Second-order centered finite differencing is used for the spatial derivatives. We introduce two physics-based preconditioners associated with the first- and second-order temporal discretizations. The theoretical time and spatial accuracies of the numerical scheme are verified through convergence tables and graphs that correspond to different computational settings. We present efficiency studies with and without using the preconditioners. Our numerical findings indicate the excellent performance of the newly proposed preconditioning strategies. In other words, when we turn the preconditioners on, the average number of GMRES and the Newton iterations are significantly reduced.

A Modified JFNK for Solving the HTR Steady State Secondary Circuit Problem

A nuclear power plant is a complex coupling system, which features multi-physics coupling between reactor physics and thermal-hydraulics in the reactor core, as well as the multi-circuit coupling between the primary circuit and the secondary circuit by the shared steam generator (SG). Especially in the pebble-bed modular HTR nuclear power plant, different nuclear steam supply modules are further coupled together through the shared main steam pipes and the related equipment in the secondary circuit, since the special configuration of multiple reactor modules connects to a steam turbine. The JFNK (Jacobian-Free Newton–Krylov) method provides a promising coupling framework to solve the whole HTR nuclear power plant problem, due to its excellent convergence rate and strong robustness. In this work, the JFNK method was modified and applied to the steady-state calculation of the HTR secondary circuit, which plays an important role in simultaneous solutions for the whole HTR nuclear power plant. The main components in the secondary circuit included SG, steam turbine, condenser, feed pump, high/low-pressure heat exchanger, deaerator, as well as the extraction steam from the steam turbine. The results showed that the JFNK method can effectively solve the steady state issue of the HTR secondary circuit. Moreover, the JFNK method could converge well within a wide range of initial values, indicating its strong robustness.

Jacobian-free Newton–Krylov method for the simulation of non-thermal plasma discharges with high-order time integration and physics-based preconditioning

A preconditioning framework for the numerical simulation of non-thermal streamer discharges is developed using the Jacobian-free Newton-Krylov (JFNK) method. A reduced plasma fluid model is considered, consisting of electrons, one positive ion, one negative ion, and the electrostatic potential. The plasma kinetics model includes ionization, electron-ion recombination, electron attachment, electron detachment, and ion-ion recombination. The governing equations are made dimensionless, discretized in space with finite differences, and integrated in time with a fully implicit method based on high-order backward differentiation formulas. The preconditioning framework is based on a linearized form of the governing equations and physics-based operator splitting. The efficiency of the preconditioning strategy is assessed through two test cases: streamer propagation between parallel plates and an axisymmetric pin-to-pin discharge. The fully implicit approach overcomes traditional restrictions in the time step size due to processes such as electron drift, electron diffusion, and dielectric relaxation. Excellent performance is observed through relevant statistics of the JFNK solver, although the number of linear iterations increases for the pin-to-pin discharge when nonlinear numerical boundary conditions are imposed at the electrodes. Performance studies show scalability with O(100-1000) processors for O(10M) unknowns with ample room for optimization.

Operator split, Picard iteration and JFNK methods based on nonlinear CMFD for transient full core models in the coupling multiphysics environment

Operator Split (OS), Picard iteration and Jacobian-free Newton Krylov (JFNK) methods based on nonlinear coarse mesh finite difference (CMFD) with nodal expansion methods (NEM) are respectively introduced to solve the transient full core coupled models in the COupling Multiphysics Environment (COME). In the COME, in addition to the traditional OS and Picard coupled methods, we develop two different JFNK-based coupled methods with physics and algebraic-based hybrid preconditioners at each time step: OS_JFNK and JFNK coupled methods. Specifically, the efficient JFNK method is applied to solve CMFD discrete models of only neutronics in the OS framework for OS_JFNK coupled strategies or whole Neutronics-Thermal Hydraulic (N-TH) coupled models for JFNK coupled methods. Then nonlinear corrective coupling coefficients are updated by the two-node NEM every few Newton steps at each time step. Six cases of NEACRP 3D core transient benchmarks are analyzed and numerical solutions of OS, OS_JFNK, Picard and JFNK methods in the COME generally agree well with other published solutions. JFNK and Picard coupled methods can track the reference solutions better than OS and OS_JFNK coupled methods due to multiple iterations or simultaneous solutions between neutronics and TH models in the tight coupled form. Furthermore, OS_JFNK or JFNK methods can obtain higher computational efficiency than OS or Picard methods with or without the fission source extrapolation acceleration, respectively, which indicates the potential of the developed COME using different coupled methods for the transient full core multiphysics coupled problems.

A modified JFNK with line search method for solving k-eigenvalue neutronics problems with thermal-hydraulics feedback

The k-eigenvalue neutronics/thermal-hydraulics coupling calculation is a key issue for reactor design and analysis. Jacobian-free Newton-Krylov (JFNK) method, featured with super-linear convergence rate and high efficiency, has been attracting more and more attention to solve the multi-physics coupling problem. However, it may converge to the high-order eigenmode because of the multiple solutions nature of the k-eigenvalue form of multi-physics coupling issue. Based on our previous work, a modified JFNK with a line search method is proposed in this work, which can find the fundamental eigenmode together with thermal-hydraulics feedback in a wide range of initial values. In detail, the existing modified JFNK method is combined with the line search strategy, so that the intermediate iterative solution can avoid a sudden divergence and be adjusted into a convergence basin smoothly. Two simplified 2-D homogeneous reactor models, a PWR model, and an HTR model, are utilized to evaluate the performance of the newly proposed JFNK method. The results show that the performance of this proposed JFNK is more robust than the existing JFNK-based methods.

Jacobian-free Newton Krylov coarse mesh finite difference algorithm based on high-order nodal expansion method for three-dimensional nuclear reactor pin-by-pin multiphysics coupled models

The nuclear reactor core is a large-scale, complicated and tight-coupled nonlinear multiphysics coupled system including neutron transport, heat conduction of nuclear fuel pins, coolant flow, heat transfer and so on. To meet the demands of the advanced nuclear reactor analysis and design, there is an urgent need to reduce the computational costs and improve the convergence rate for the three-dimensional (3D) nuclear reactor core multiphysics coupled models, especially for a high-fidelity pin-by-pin (at the level of fuel pin by fuel pin) full core simulation. In this work, we develop CNJ, an efficient Jacobian free Newton Krylov (JFNK) coarse mesh finite difference (CMFD) algorithm based on the nodal expansion method (NEM) to simultaneously solve the complicated pin-by-pin neutronics-thermal hydraulic (N-TH) coupled models. The high-order NEM can improve the numerical accuracy of the CMFD method on a coarse mesh and the JFNK method can ensure the rapid convergence and high efficiency on large-scale and nonlinear coupled discrete systems. By combining the CMFD, NEM and JFNK methods and making full use of their respective advantages, the CNJ algorithm can obtain high accuracy and efficiency even on a coarse mesh to greatly reduce the number of solution variables and the computational cost of the 3D coupled models. Finally, two representative and complicated 3D pin-by-pin N-TH coupled models are analyzed to evaluate the numerical accuracy and efficiency advantage of CNJ in detailed.

On the Jacobian approximation in sea ice models with viscous-plastic rheology

Viscous-plastic rheology is widely used in the sea ice modeling community at increasingly high resolutions. Due to the high degree of non-linearity of the rheological constraints, accurate approximation of the Jacobian is required to improve the efficiency of the implicit solvers of the sea ice momentum equation in pack ice. We consider the analytical Jacobian of the ice momentum equation and assess its approximation errors in the Jacobian-free Newton–Krylov (JFNK) method and in the family of more traditional schemes which neglect the dependence of viscosity coefficients on the deformation tensor. It is shown that this dependence provides a substantial contribution to the Jacobian, especially in the regions enriched by linear kinematic features like ridges and elongated polynyas. Numerical experiments indicate that performance of the Newton solvers is also sensitive to errors associated with inexact computation of the search direction, that may be caused, in particular, by numerical approximations of the Jacobian which violate its dissipative property. Based on this analysis, an improved selective damping strategy for the Newtonian solver of the momentum equation is proposed. A series of numerical experiments conducted in simulated pack ice environment demonstrate faster convergence of the updated solver with analytical Jacobian as compared to the one based on the selectively damped JFNK method with inexact GMRES solver in the inner loop.

An effective method for modeling 1D two-phase two-fluid six-equation model with automatic differentiation approach

The modeling of the two-fluid model with plenty of realistic closure correlations has been a tough problem because of the difficulty of forming the Jacobian matrix manually. To predict the phase appearance/disappearance in two-phase flow problems, the high-order spatial and fully implicit temporal discretization typically should be adopted. However, considering the requirement of a fairly accurate evaluation of the Jacobian matrix, the implementation of high-order schemes is a challenging task. In this study, the automatic differentiation (AD) method was employed to overcome these difficulties, which seems the first demonstration of the high-order scheme and the AD method in the two-fluid model with realistic relations. With the AD method, the high-order fully implicit scheme, and the preconditioned Jacobian-free Newton-Krylov (PJFNK) method, the two-fluid code was developed to solve the two-phase flow and heat transfer problems in this study. Through simulations of a series of hydrodynamics two-phase flow problems, this code was verified and validated in both transients and boiling conditions. Meanwhile, code-to-code comparisons with well-validated code, RELAP5, were presented. The results demonstrate the promising performance of the developed code for solving 1D two-phase flow problems.

A multi-level nonlinear elimination-based JFNK method for multi-scale multi-physics coupling problem in pebble-bed HTR

In the pebble-bed high-temperature gas-cooled reactor, various physical phenomena occur at different spatial scales and are coupled together tightly in a nonlinear manner. More specifically, the local temperature profiles inside the fuel sphere and TRISO particle, which are essential in reactor design and safety analysis, are coupled with the global physical fields, such as the pebble bed temperature, helium temperature, as well as the fission heat source. Therefore, the multi-physics calculation combined with multi-scale model should be employed to describe the reactor behavior. The computational efficiency is the critical issue for the numerical solution of such global and local coupled systems. Jacobian-free Newton Krylov (JFNK) algorithm is an excellent candidate for multi-physics coupling issues due to its high convergence rate and superior computational performance, but it doesn’t consider the physical features at different scales. In this work, a multi-level nonlinear elimination is proposed within the JFNK framework to deal with the multi-scale coupled system in high-temperature gas-cooled reactor. The numerous local multi-scale temperature variables, including the fuel sphere and TRISO particle temperatures, are eliminated efficiently based on the Schur complement factorization due to its well sparse matrix structure. As a result, the number of unknowns is significantly reduced and the computational performance is improved. A steady-state pebble-bed HTR model is used in this work, and the numerical results show that multi-level nonlinear elimination-based JFNK method has high computational efficiency and numerical robustness, whose performance is 4 ∼ 6 times higher than that of the original JFNK algorithm.

Large-Eddy Simulation of Transonic Buffet Using Matrix-Free Discontinuous Galerkin Method

We present an implicit large-eddy simulation of transonic buffet over the OAT15A supercritical airfoil at Mach number 0.73, angle of attack 3.5 deg, and Reynolds number . The simulation is performed using a matrix-free discontinuous Galerkin (DG) method and a diagonally implicit Runge–Kutta scheme on graphics processor units. We propose a Jacobian-free Newton–Krylov method to solve nonlinear systems arising from the discretization of the Navier–Stokes equations. The method successfully predicts the buffet onset, the buffet frequency, and turbulence statistics owing to the high-order DG discretization and an efficient mesh refinement for the laminar and turbulent boundary layers. A number of physical phenomena present in the experiment are captured in our simulation, including periodical low-frequency oscillations of shock wave in the streamwise direction, strong shear layer detached from the shock wave due to shock-wave/boundary-layer interaction and small-scale structures broken down by the shear-layer instability in the transition region, and shock-induced flow separation. The pressure coefficient, the root mean square of the fluctuating pressure, and the streamwise range of the shock wave oscillation agree well with the experimental data. The results suggest that the proposed method can accurately predict the onset of turbulence and buffet phenomena at high Reynolds numbers without a subgrid scale model or a wall model.

An adaptive scalable fully implicit algorithm based on stabilized finite element for reduced visco-resistive MHD

The magnetohydrodynamics (MHD) equations are continuum models used in the study of a wide range of plasma physics systems, including the evolution of complex plasma dynamics in tokamak disruptions. However, efficient numerical solution methods for MHD are extremely challenging due to disparate time and length scales, strong hyperbolic phenomena, and nonlinearity. Therefore the development of scalable, implicit MHD algorithms and high-resolution adaptive mesh refinement strategies is of considerable importance. In this work, we develop a high-order stabilized finite-element algorithm for the reduced visco-resistive MHD equations based on the MFEM finite element library (mfem.org). The scheme is fully implicit, solved with the Jacobian-free Newton-Krylov (JFNK) method with a physics-based preconditioning strategy. Our preconditioning strategy is a generalization of the physics-based preconditioning methods in Chacón et al. (2002) [3] to adaptive, stabilized finite elements. Algebraic multigrid methods are used to invert sub-block operators to achieve scalability. A parallel adaptive mesh refinement scheme with dynamic load-balancing is implemented to efficiently resolve the multi-scale spatial features of the system. Our implementation uses the MFEM framework, which provides arbitrary-order polynomials and flexible adaptive conforming and non-conforming meshes capabilities. Results demonstrate the accuracy, efficiency, and scalability of the implicit scheme in the presence of large scale disparity. The potential of the AMR approach is demonstrated on an island coalescence problem in the high Lundquist-number regime (≥107) with the successful resolution of plasmoid instabilities and thin current sheets.