High Performance Computing for MHD Ramakanth Munipalli, P.-Y. Huang HyPerComp Inc., Westlake Village, CA 91361 In collaboration with S.Smolentsev, N.B.Morley, A.Ying, G.Pulugundla, M.Abdou, UCLA 2nd EU-US DCLL Workshop, UCLA, November 15, 2014 Research supported by current and prior SBIR funding from DOE Outline 1. HyPerComp Inc. – who are we? Small company, located 30 miles NW of LA, estd. 1998, Lengthy, fruitful collaboration with UCLA in numerical MHD 2. MHD – what have we done so far in fusion related MHD modeling? 3. Multiphysical extensions Mass transport, corrosion, integrated multiphysical modeling 4. What are we working on now? 5. High performance computing – some tools and techniques Some prospects and future directions 1. HyPerComp Inc. – Who are we? HyPerComp develops high performance computing software in fluid mechanics and electromagnetics 3 Recent Examples: Very large scale CPU/GPU parallel electromagnetics modeling Liquid rocket combustion dynamics Coupled fluid-structure-acoustics modeling for rotorcraft A Technology Portfolio Computational Electromagnetics RCS, SAR Imagery and range profiles, Penetrable materials, Antennas, Optics Magnetohydrodynamics & Plasmadynamics Flow Control, Liquid metal MHD – nuclear fusion, Combustion control Computational Fluid Dynamics Incompressible to Hypersonic flow, Turbulent Combustion, Rotorcraft, Acoustics, Morphing structures, Optimization High Performance Multiphysical Computations Very high order accurate solutions, CPU/GPU based solvers, Reduced Basis Models TEMPUS HOME HD Physics HYCE 2. MHD – what have we done so far? MHD Top-level concern in fusion blanket studies (pressure drop, corrosion, tritium transport, etc.) Difficult problem to solve No analytical solutions in all but the most elementary flows Numerical solutions riddled with difficulties and uncertainties HyPerComp Incompressible MHD solver for Arbitrary Geometry HIMAG was developed by HyPerComp in a research partnership with UCLA to model the flow of liquid metals in nuclear fusion reactor design. These flows are characterized by very high magnetic field strength, complex geometry and fluid-solid coupling via electric current and heat. HIMAG is a pioneer in the reliable computation of such flows among both commercial as well as research simulation software. HIMAG: Technical Summary HIMAG is a parallel, unstructured mesh based MHD solver. High accuracy at high Hartmann numbers is maintained even on nonorthogonal meshes HIMAG can model single-phase as well as two-phase (free surface) flows Multiple conducting solid walls may be present in the computational domain Graphical User Interfaces are provided for the full execution of HIMAG Heat transfer, natural convection, temperature dependent properties can be modeled Extensive validation and benchmarking has been performed for canonical problems. Cases involving Ha > 1000 have never been demonstrated on non-rectangular meshes prior to HIMAG Governing Equations ui xi t Two phase fluid flow without MHD 0, u i xi 0, 2 1 2 H ( ) ui t u i u j xj 2 1 2 H ( ) , ui xi x j x j p Two phase fluid flow with MHD, heat transfer ui t u i u j t 1 u i J B xi x i x i 1 p viscous T t ui T xi 1 cp k T Lorentz J E Q g gravity x surface tension T T ref g buoyancy Treatment of electromagnetic fields in liquid metal MHD When the magnetic Reynolds number is low (Rem = μ0σLU << 1), induced magnetic field is negligible. Applied magnetic field is the total magnetic field, and an electric potential may be derived: Ohm’s law: J V B E The potential is obtained from the condition that the current density is solenoidal J 0 V B A more general treatment solves for the induced magnetic field: B t V B B V 1 2 B in the fluid, J B Electric current can flow across interfaces between media solid fluid Appropriate simplifications (J=0, etc.), or, Boundary conditions (normal current = 0) can be used Numerical difficulties in computing high Ha MHD flows 1. The need to resolve and minimize numerical errors in Hartmann layers (Thickness about 1/Ha) 2. The pace of convergence of Poisson equation solvers (for Pressure, Electric potential, divergence of B) 3. Computationally intensive corrections for non-orthogonal meshes 4. Long periods of integration needed to account for flow development, unsteady effects 5. Time taken to develop CAD to CFD mesh for high Ha problems Flow enters a strong B-field The high Hartmann number problem As Ha = BL increases, regions of sharp variations in velocity and electric current density appear in the flow. Numerically, small inaccuracies in the Hartmann layer are greatly amplified The exacting needs of numerical MHD FLOW B S id e la ye rs ~ 1 Ha 1 / sqrt(Ha) H a rtm a n n la ye rs ~ B 1 U(z) Ha j(y,z) 2h d//≈Ha-1/2 Side layer 1 / Ha dHa≈Ha-1 Hartmann layer Verification against canonical benchmark problems Fully developed flow at Ha = 10,000 in a duct with square cross section – compared with the exact analytical solution Fully developed flow at Ha = 1,000 in a duct with circular cross section Natural Convection A good match of Nusselt number with published data has been observed. Second order accuracy has been ascertained. Ra 103 104 105 106 Ra = 103 Mesh NuHIMAG Nucomp 21x21 41x41 81x81 1.1227 1.1191 1.1181 1.118 21x21 41x41 81x81 2.3084 2.2611 2.2488 2.245 21x21 41x41 81x81 4.9370 4.6312 4.5490 4.522 21x21 41x41 81x81 161x161 10.661 9.4598 8.9885 8.8659 Ra = 104 Ra = 105 Ra = 106 8.829 Validation against experimental data: ALEX duct experiment (1987, ANL) Circular duct Ha=5800 Square duct HIMAG: Some MHD specific numerics Contact resistance Wall functions z dH ~ S olid-fluid Interface w ith contact resistance 1 Ha u f u ic 1 d Ha e N Ha N B f s =0 = 10 1 2 y j s f K J n Coarse mesh in Ha layers Jump in potential at arbitrary material interfaces can now be captured Robust modeling of strong natural convection Newton-Krylov based schemes are used to perform matrix inversion in a upwind semi-implicit procedure to stabilize the simulation of flow with strong natural convection (Gr = 10 9). Flow in a 3-D channel with MHD and heat transfer (results in next chart) Gr = 109 in a square cavity Streamlines (left) and isotherms (right) X B 2 g Ym 2a=20 cm 40 cm A case study in the use of wall functions in solution acceleration Flow Re = 10,000, Ha = 400, Gr = 107 Wall functions are used to model MHD flow in an insulating channel with dimensions similar to DCLL. These functions are applied at all walls (iwall=10) or only at Hartmann walls (iwall=11) Full MHD solution (left) and solution with wall functions (right) Sample speedup results on 16 CPUs Upward flow Full solution iwall=0 g Temperature # cells Speedup 299,440 1 Wall function iwall=10 70,080 19.18 Wall function iwall=11 162,336 6.2 Velocity contours Comparison of full and wall function solutions along centerline DCLL concept and crucial MHD features T o p P la te FW B a ck P la te A ssy S iC G rid P la te A ssy B o tto m P la te Full 3-D MHD (with natural convection) model of the DCLL B z y x g Flow enters from lower right, exits above A low conductivity flow channel insert is present Velocity profiles shown on right and pressure on left 3. Multiphysical studies Fusion relevant liquid metal MHD includes a variety of multiphysical phenomena Fluid, heat and mass transport Natural convection Steady and unsteady electromagnetic phenomena Contact resistance (thermal and electric) Ferromagnetic effects Two phase flow, surface tension Phase change (due to high heat flux, pressure variation, etc.) Corrosion, electrochemistry at walls Fluid-structure interaction HIMAG: Free surface flow simulations MTOR, experiment at UCLA Liquid metal jet in a magnetic field HIMAG: Free surface flow simulations – contd. Plasma-liquid metal interaction in a magnetic field: DiMES Plasma Plasma-liquid interaction Liquid metal HIMAG uses the level set method in unstructured meshes. Method permits large deformations of the free surface. We are presently working on massively parallel free surface capture simulations with scalable adaptive meshing (w. RPI) Mass Transport Tritium transport We built an independent software named CATRIS from the basic data structure of HIMAG, to focus on specific issues in mass transport. Corrosion Some capabilities to simulate corrosion were added to CATRIS, together with Lagrangian models for particulate transport CATRIS (Corrosion And TRItium transport Solver) Written as a new stand-alone system, CATRIS focuses on the following: (a) the transport of tritium and its permeation through walls (b) corrosion and deposition of iron contained within structural materials of the system. Blanket Region Individual components from outer region Whole system modelling CFD/MHD Solver Thermal Hydraulics Code Velocity and Temperature CATRIS Code Tritium Transport Preprocessing Grid Generation Dissolution and diffusion Properties Transport in solids and fluids Interfacial phenomena Effect of He bubbles Corrosion Dissolution into PbLi Surface layer effects Nucleation Transport of Fe particles Deposition in cold regions Numerical solvers Advection/diffusion equation Two-fluid model Particle tracking Visualization Tritium transport and particle transport models in CATRIS DCLL Geometry (not to scale) Tritium concentration in module computed for B=4T (Ha = 15 000), u=0.065m/s. B z y x RAFS wall 5 mm thick 2 mm gap y 2.0 m 2.26 m Outflow z 211 mm FCI SiC wall 5 mm thick Inflow 207 mm 0.3 m 231 mm A sample study to inject particles in an MHD flow, which migrate in response to the gradient in an applied magnetic field. These particles will be produced at walls using a corrosion BC and deposited likewise. 4. What are we working on now? Current Research Objectives Brand new software implementation of the induced magnetic field formulation – adds ability to model high magnetic Reynolds number, strongly unsteady flows Improvement in time accuracy Enhanced robustness and speed – better parallelization Dramatically reduce simulation time for MHD flows in geometries of practical interest – DCLL, etc. Transition to general EM applications in materials processing Induced magnetic field formulation 1 B u B t m B • Boundary conditions for the induced magnetic field b: Approximate boundary condition: b B B0 b || b b n n 0 b n 0 n Magnetic domain boundary condition b 0 b 0 n Functional Decay based boundary conditions 1 b b r r – distance Integro-different boundary conditions b so b 0 b\\ bn Case 1 – 3-D Lid-driven cavity with conducting walls except top lid (Re=100, Ha=10, Cw=0.4) U B a) Schematic view of 3-D lid-driven cavity with conducting walls (Blue area), Red – Liquid. b) Velocity profile along x=0 and z=0 sections Case 1 – Cavity with conducting walls - Continued e) Comparison of velocity curves by B-formulation with those by potential solver (HIMAG). f) Convergence history. Case 2 – 3-D Lid-driven cavity with insulating walls (Re=100, Ha=45) U B a) Schematic view of 3-D liddriven cavity with insulating walls: Red – Liquid. b) Velocity profile along x=0.5 and y=0.5 sections Case 2 – Cavity with insulating walls - Continued Comparison of velocity w(x) by Bformulation with published data Comparison of velocity u(z) by Bformulation with published data 5. High Performance Computing Summary of our approach Mathematical methods Multigrid methods: Agglomeration, unnested premeshed scheme, algebraic multigrid Hybrid meshes Implicit schemes – time stepping, accuracy CN, AB, SIRK schemes SIMPLE scheme for steady flow Local time stepping Full matrix solvers – BICGSTAB on CPU/GPU Interpolation procedures Data storage: Non-orthogonal correction GPU-based programming Rapid prototyping Canonical decomposition of DCLL geometry User interfaces for inlet, fci, bend, etc. Template for multigrid (3-D) – dir. agglomeration Template for hybrid mesh generation blocking Applications of analysis EM coupling between neighboring channels Fringing fields, self-consistent field formulations Wall functions – insulating, perfectly conducting Wall functions for fringing fields Patched analytical – numerical solution Combination of duct flow solutions Poisson solver acceleration by sparse matrix storage and inversion Cases Poisson solver with Neumann BCs option Run time Log L2 norm CG iortho=1 95.7 s -0.479 CG, K saved iortho=1 57.9 s -0.479 CG iortho=2 221.9 s -1.958 CG, K saved iortho=2 138.9 s -1.963 Cases Hunt’s fully developed channel flow option Run time CG iortho=1, nmhd=1 298 s CG, K saved iortho=1, nmhd=5 168 s Option Run time iortho=2, nppe=1, nmhd=1 29783 s iortho=2, nppe=5, nmhd=5 19633 s Option Run time CG iortho=1, nppe=1 19996 s CG, K saved iortho=1, nppe=5 13060 s Cases option Run time CG iortho=1, nppe=1, nmhd=1 3205 s CG, K saved iortho=1, nppe=5, nmhd=5 1712 s Cases 3-D circular duct Ha = 10, CG Re = 10, cw = 0.1 CG, K saved Cases 3-D rectangular duct Ha = 0, Re = 15,250 3-D rectangular duct Ha = 100, Re = 10, cw = 0.1 The use of non-orthogonal meshes Among the many benefits of the HIMAG approach has been the preservation of accuracy on arbitrary mesh systems – making the calculation more efficient Hybrid/unstructured meshes reduce number of mesh points needed, by transitioning between regions of different resolution requirements “Template” driven mesh generation HIMAG: A nested multigrid method A sequence of unnested grids for multigrid Poisson solution, showing reduced convergence time (above), automatic generation of a sequence of meshes from CAD input (below) Some prospects and future directions Summary We have a number of ongoing software development activities which can serve simulation needs in fusion: Robust, comprehensive, self-consistent MHD physics modeling CPU/GPU parallel computing Highly scalable parallel mesh adaptation Strengths in coupled MHD/heat/mass transfer analysis (numerous existing physical and numerical models) Integrated and customizable software solutions Future Directions EM toolkit for flow and materials processing 1. Combined calculation of fluid flow, heat transfer (radiative, convective), mass transfer (including corrosion) and electromagnetics - applications in aerospace propulsion, materials processing 2. Localized plasma models – thermal as well as weak ionization, extended to large length scales as source terms (external flow control, plasma assisted processing, radar analysis) 3. Two phase flow - models of free surfaces such as melt layers 4. Phase change - evaporation, melting, solidification, crystallization, etc.