Mixed Generalized Multiscale Finite Element Method for a Simplified Magnetohydrodynamics Problem in Perforated Domains

In this paper, we consider a coupled system of equations that describes simplified magnetohydrodynamics (MHD) problem in perforated domains. We construct a fine grid that resolves the perforations on the grid level in order to use a traditional approximation. For the solution on the fine grid, we construct approximation using the mixed finite element method. To reduce the size of the fine grid system, we will develop a Mixed Generalized Multiscale Finite Element Method (Mixed GMsFEM). The method differs from existing approaches and requires some modifications to represent the flow and magnetic fields. Numerical results are presented for a two-dimensional model problem in perforated domains. This model problem is a special case for the general 3D problem. We study the influence of the number of multiscale basis functions on the accuracy of the method and show that the proposed method provides a good accuracy with few basis functions.


Introduction
The magnetohydrodynamic (MHD) system describes the interactions of electrically conducting incompressible flows in the presence of a magnetic field. MHD flows are usually involved in the contexts: liquid metal cooling of nuclear reactors, electromagnetic casting of metals, MHD generators, accelerators, and MHD ion propulsion, see [1][2][3]. The equations that model the MHD are the coupling of incompressible Navier-Stokes equations and the Maxwell's equations via the Lorentz force and Ohm's law. There are some research works devoted to the numerical approximation of MHD equations. For treating the nonlinear terms effectively, three classical iterative methods (Stokes-type, Newtown and Oseen-type ones) within finite element approximation for the steady MHD equations have been developed in [4,5]. Especially, the theoretical analysis and numerical examples in [4,6] show that the Picard iteration (called Oseen-type iteration by the authors in the reference) is suitable for high Reynolds number problems. The local and parallel finite element algorithms based on two-grid methods were given in [7][8][9]. The stabilized finite element schemes that capture physical solutions were investigated in [10]. The mixed finite element methods with divergence-free velocities and magnetics fields are developed [11][12][13]. However, analysis of the above literature mainly focuses on the singly-connected computational domains. In this paper, our goal is to develop multiscale methods for MHD flows on complex domains with perforations.
In this paper, we consider a simplified MHD problem in perforated domains, and propose a multiscale method for constructing the reduced order model on the coarse grid. In particular, we consider a special 2D case of MHD equations and the full MHD problem will be considered in the future. We start with the mixed formulation for both equations of the simplified MHD problem in a perforated domain. For the solution of the nonlinear system, we use a Picard iteration to decouple the equations. We construct a fine grid that resolves perforations on the grid level and write an approximation using a mixed finite element method. Problems in perforated domains have a multiscale nature and approximations on the fine grid lead the large system of equations. To reduce the size of the discrete system, multiscale methods or homogenization techniques are needed [14][15][16][17][18][19]. Homogenization techniques are constructed to solve problems with scale separation, i.e., when perforations are of similar sizes. When the perforated media do not have scale separation, multiscale methods are used [20][21][22][23][24]. In the papers [20,21], the authors presented the multiscale finite element method for problems in perforated domains, but these approaches use a limited number of degrees of freedom per coarse element. In the papers [18,25], the general approach for constructing the multiscale basis functions using local spectral characteristics was presented in [18,25]. Continuous and Discontinuous Galerkin Generalized Multiscale Finite Element Methods (GMsFEM) for coarse grid approximation of the problems in perforated domains was presented in [26,27]. The GMsFEM is a systematic approach to identify multiscale basis functions via local snapshots and local spectral problems. The local snapshots are constructed by solving local problems and contain the information about local geometry structure (perforations). By performing local spectral decomposition, the method identifies multiscale basis functions. For accurately approximating the fluxes, we use a Mixed Generalized Multiscale Finite Element method and construct multiscale basis functions for them [28][29][30][31]. In Mixed GMsFEM, we construct the snapshot space for each coarse edge, perform the local spectral decomposition in the snapshot space and select the dominant modes as multiscale basis functions for the flux. The presented work is based on our previous papers [30,31], where we presented a Mixed GMsFEM for solutions of the problems in perforated domains. In this paper, we apply and study the method for solving the incompressible magnetohydrodynamics (MHD) problem in perforated domains. We present numerical results for some perforated geometries 2D. These numerical results demonstrate one can achieve a good accuracy with fewer basis functions.
The paper is organized as follows. In Section 2, we describe the problem formulation. Next, we decouple the system of equations using Picard iterations. In Sections 3 and 4, we present a fine grid and multiscale approximations for the magnetic field and flow problems, respectively. Fine grid and coarse grid approximations are presented using mixed formulation. Construction of the multiscale basis functions are presented. In Section 5, we present numerical results for the two-dimensional problem in a perforated domain. The paper ends with a conclusion.

Problem Formulation
The full steady incompressible MHD equations in Ω ⊂ R 3 is a complex system which link the velocity field u = (u 1 , u 2 , u 3 ), magnetic field B = (B 1 , B 2 , B 3 ) and pressure p: where R e is the hydrodynamic Reynolds number, R e m the magnetic Reynolds number, S c the coupling number. We refer the readers to [32] for details on the existence, uniqueness, and finite element approximation of solutions of (1).
This article considers a magnetic field B = (0, 0, B) perpendicular to the velocity field u = (u 1 , u 2 , 0) of (1) in two-dimensional perforated domain Ω, then the simplified MHD equations without any body force and hydrodynamic convective terms are described as follows: Let D = R −1 e m , q = −D∇B be the flux and Ω is the perforated domain (see Figure 1). We write Equation (2) in the mixed form and consider with the following boundary conditions where Γ P is the boundary of perforations, Γ 1 ∪ Γ 2 ∪ Γ P = ∂Ω, n is the unit outward normal vector on ∂Ω and I is the n × n identity matrix. For linearization, we use Picard iterations and obtain the following algorithm Here k is the nonlinear iteration.

Magnetic Field Problem
We have the following variational formulations: where k is the nonlinear iteration and We can rewrite fine grid approximation in the matrix form where C k = C(u k ). The fine grid problem is solved using the lowest-order Raviart-Thomas element.

Multiscale Approximation
To construct the multiscale space for the coarse scale approximation, we use a Mixed GMsFEM [23,[28][29][30][31]. Let T H be a coarse grid of domain Ω and E H be the set of all faces of the coarse grid, where H is the coarse mesh size and N e be the total number of faces. We define the neighborhood of the face E i ∈ E H by and ω i is a union of two coarse grid blocks when E i lies in the interior of the domain Ω (see Figure 2). We construct a multiscale space for the flux and for the B, we use the space of piecewise constant functions over the coarse triangulation T H .
Here ψ i are the multiscale basis functions that supported in a coarse neighborhood For the coarse grid approximation, we use a Mixed GMsFEM and have the following variational formulation: We can write our system in matrix form for each nonlinear iteration where Here R q and R B are the projection matrices We remark that {ψ i } For construction of the multiscale space for the flux, we start with the construction of the snapshot space that contains an extensive set of basis functions formed by the solution of local problems with all possible boundary conditions up to the fine grid resolution. After that, we solve a spectral problem to select dominant modes in the snapshot space.
Snapshot space. We solve the following problem on the coarse neighborhood ω i that corresponds to the coarse-grid edge with boundary condition φ i l · n i = 0 on ∂ω i and additional boundary condition φ i l · n i = δ i l on coarse edge E i , where n i is the outward unit-normal vector on ∂ω i . The local problem (10) is solved separately on each coarse-grid element K i j ⊂ ω i using the fine-grid defined in ω i by the lowest-order Raviart-Thomas element. Here E i = ∪ J i l=1 e l , where e l are the fine-grid edges on E i /Γ P , and J i is the number of the fine grid edges on E i /Γ P . The delta function δ i l is a piecewise constant function defined on E i /Γ P such that it has value 1 on e l and value 0 on the other fine-grid edges. The constant c in (10) is chosen by compatibility condition, c = 1 The collection of the solutions of the above local problems generates the snapshot space in ω i Multiscale space. For reduction on the snapshot space, we solve the following local spectral where and To form a multiscale space, we arrange the eigenvalues in increasing order, and choose the first M i eigenvalues and take the corresponding eigenvectors ψ i k = R i,snap ψ i,snap k as the basis functions (k = 1, ..., M i ) This space will be used as the approximation space for the flux.

Flow Problem
For the approximation of the flow problem on the fine grid, we use a Discontinuous Galerkin method [33][34][35][36][37][38]. Let T h be a fine-grid partition of the domain Ω, E h be the set of facets in T h and E h = E h int ∪ E h out (E h int is the set of interior facets and E h out is the set of boundary facets). We use the notations K and E to denote cell and facet in the fine grid T h . We define the jump [u] and the average {u} of a function u on interior facet where n is the unit normal vector on E, K + and K − are the two cells sharing the facet E.
We have the following variational formulations for the DG approach: find (u k+1 , p n+1 ) ∈ V u × Q p such that where The fine scale velocity space V u = {v ∈ L 2 (Ω) : v| K ∈ (P 1 (T)) 2 , ∀K ∈ T h } contains functions that are piecewise linear in each fine-grid element K and are continuous along the fine-grid edges, but are discontinuous across coarse grid edges. For the pressure, we use the space of piecewise constant functions.
In the matrix form, we have where F k+1 u = F u (B k+1 , q k+1 ).
Next, we present a multiscale method for solutions of the Stokes problem. Similarly, to the Mixed GMsFEM to obtain multiscale basis functions, we solve problems in local domains with various boundary conditions to form a snapshot space and use a spectral problem to perform a dimension reduction.

Multiscale Approximation
Next, we present the construction of the multiscale space for the coarse scale approximation of velocity. Note that, in the DG multiscale method we construct a multiscale basis in each coarse block [26,27,[39][40][41]. Let T H be a coarse-grid partition of the domain Ω with mesh size H and E H be the set of facets in T H , E H = E H int ∪ E H out . For the pressure approximation, we use the piecewise constant function space Q p,H over the coarse cells. We define V u,H as the multiscale velocity space, which contains a set of basis functions supported in each coarse block K (see Figure 3). We construct a multiscale space for the velocity T H V u,H = span{ψ i } N u i=1 , Q p,H = {r ∈ L 2 (Ω) : r| K ∈ P 0 (K), ∀K ∈ T H } and for the pressure, we use the space of piecewise constant functions over the coarse triangulation, N u = dim(V u,H ) is the number of basis functions and N p = dim(Q p,H ) is equal to the number of coarse grid cells. For the coarse grid approximation, we use a DG approach and have the following variational formulation: find (u k+1 We can write the system in matrix form for each nonlinear iteration where Here R u and R p are the projection matrices are the multiscale basis functions for velocity. For construction of the multiscale space for the velocity, we start with the construction of the snapshot space that contains an extensive set of basis functions formed by the solution of local problems with all possible boundary conditions up to the fine grid resolution. After that, we solve a spectral problem to select dominant modes of the snapshot space. Snapshot space. We construct local snapshot basis in each coarse block K i , (i = 1, · · · , N), where N is the number of coarse blocks in Ω. The local snapshot space consists of functions which are where J i is the number of fine grid nodes on the boundary of K i , and δ l i is the discrete delta function defined on ∂K i . This problem is solved on the fine mesh using some appropriate approximation spaces (DG). Here constant c is chosen by the compatibility condition,

Multiscale space.
To reduce the size of the snapshot space, we solve the following local spectral problem in the snapshot space for K i where and A i is the matrix representation of the bilinear form a i (u, v) and S i is the matrix representation of the bilinear form s i (u, v) We remark that the integral in s i (u, v) is defined on the boundary of the coarse block. In this case, the number of the spectral problem equals the number of coarse blocks. We arrange the eigenvalues in increasing order and choose the first eigenvectors corresponding to the first M i the smallest eigenvalues ψ i k = R i,snap ψ i,snap k as the basis functions (k = 1, ..., M i ) This space will be used as the approximation space for the velocity. The global multiscale space V u,H is the combination of the local ones, i.e.,

Numerical Results
We present numerical results for the model problem in the perforated domain. We use the Generalized Multiscale Finite Element method (GMsFEM) to construct a coarse grid approximation.
The presented method based on the calculations of local multiscale basis functions in order to form a projection matrix that is used to construct a lower-dimensional approximation. After obtaining the solution of the reduced-order system, we reconstruct a fine grid resolution using a transposed projection matrix. For error calculations, we use a solution of the problem with traditional finite element approximation on the fine grid as the reference solution.
A computational perforated domain is presented in Figure 4 and has dimensions Ω = [0, 1] × [0, 1]. For the approximation of the perforated domain, we use an unstructured fine grid with triangular cells that resolves perforations on the grid level. This can lead to the large discrete system of equations for the traditional finite element approximations. To reduce the size of the system, we present a multiscale solver that used to construct a reduced-order approximation on the coarse grid. The coarse grid can be either structured or unstructured. In this work, for simplicity, we consider a 10 × 10 structured grid to illustrate the multiscale method. In many applications, it is easier to use structured coarse grids and resolve the subgrid effects as we do in our simulations. We construct an unstructured fine mesh that contains 61,912 cells, 93,668 facets and 31,727 vertices. The coarse grid is uniform, which contains 100 cells, 220 facets and 121 vertices. We set R −1 e = 1, D = 10, S c = 1 and g = 1. In Figures 5-7, we present the results of the proposed multiscale method. The distribution of the magnitude fields (B and q) are presented in Figures 5 and 6. In Figure 7, we show the velocity field. At the top of the figures, we presented the fine-scale solutions. Multiscale solutions are present on the bottom figure with DOF c = 980 for q and DOF c = 4100 for velocity field u. Note that, we use the piecewise constant functions for multiscale solution of B. We use a Picard iteration and perform three iterations for all cases with B = 10 −9 .   We calculate following errors: where (q, B) and (u, p) are the reference solutions using traditional fine grid solvers,B andp are the coarse cell average for reference solutions, (q H , B H ) and (u H , p H ) are the coarse grid solutions using presented multiscale solvers, q ms = R T q q H and u ms = R T u u H are the reconstructed multiscale solutions. In Table 1, we show the relative error for q and B. We presented errors for different velocities: (1) multiscale velocities with DOF c = 1100 (10 basis), (2) multiscale velocities with DOF c = 2100 (20 basis), (3) multiscale velocities with DOF c = 3100 (30 basis), (4) multiscale velocities with DOF c = 4100 (40 basis), (5) multiscale velocities with DOF c = 6100 (60 basis) and (6) fine-scale solution with DOF f = 433,384. We see the relative error for q is reduced from 4.1% to 0.8%, and the coarse cell average error for B is around one percent. In Table 2 we present the relative errors for velocity and for pressure field. We presented errors for u and p with different B: (1) multiscale magnetic field with DOF c = 320 (1 basis), (2) multiscale magnetic field with DOF c = 540 (2 basis), (3) multiscale magnetic field with DOF c = 760 (3 basis), (4) multiscale magnetic field with DOF c = 980 (4 basis) and (5) fine-scale solution with DOF f = 155,580. We see the relative velocity error reduces from 60% to 8% with a larger number of multiscale basis functions. We use a discontinuous Galerkin Generalized Multiscale Finite Element Method (DG-GMsFEM) for the flow problem and Mixed Generalized Multiscale Finite Element Method (Mixed GMsFEM) for the magnetic field problem. The construction of basis functions in both methods follows the general framework of GMsFEM by using local snapshots and local spectral problems. The errors can vary depending on a choice of discretization. In order to further error reduction for the multiscale approximation of the flow problem, we will apply oversampling techniques in future work [39].

Conclusions
We considered the magnetohydrodynamics (MHD) problem in a perforated domain. A mathematical model was described using a coupled system of equations for the magnetic field and for the flow problem. We constructed a fine grid that resolved perforations on the grid level in order to use a traditional approximation. For the solution on the fine grid, we constructed approximation using the mixed finite element method. To reduce the size of the fine grid system with an accurate solution on the coarse grid, a Mixed Generalized Multiscale Finite Element Method (Mixed GMsFEM) was presented. Numerical results were presented for a two-dimensional model problem in a perforated domain and we studied the influence of the number of multiscale basis functions on the method accuracy. In the multiscale solver, the main computational gain was due to fewer multiscale basis functions, which are used to reduce the size of the system. For the traditional fine grid solutions, we have DOF f = 155,580 for the magnetic field problem and DOF f = 433,384 for the flow problem. For multiscale solvers, we have 0.9% error for q with three multiscale basis functions with DOF c = 760 (0.48% of DOF f ) and 9% error for u with 40 multiscale basis functions with DOF c = 4100 (0.9% of DOF f ). In the future, we plan to consider oversampling techniques and construct a multiscale solver for three-dimensional problems.