Definite Maxwell Problem

by Samik Mukhopadhyay and Erin Okey, Brown University

45 minutes basic

Maxwell's equations are a set of coupled PDEs that form the basis for the study of electromagnetism. This annotated example demonstrates how to solve the second-order definite Maxwell equation, $$ \nabla \times (\mu^{-1} \nabla \times E) + \sigma E = f, $$ corresponding to a simple model of electromagnetic diffusion with permeability and conductivity parameters $\mu>0$ and $\sigma>0$, respectively. The boundary condition is given by setting $E \times n$ equal to some given tangential field $g$.

After taking the inner product of (1) with a test vector field $F$, followed by integrating over the domain $\Omega$, we can apply Stokes' theorem to arrive at the following variational formulation: find $E \in H(\text{curl})$ with $E \times n = g$ on the boundary of our domain $\partial \Omega$ such that $$a(E,F) := (\mu^{-1} \nabla \times E, \nabla \times F) + (\sigma E, F) = (f,F)=: b(F)$$

for all $F \in H(\text{curl},\Omega)$ with $F \times n = 0$ on $\partial \Omega$.

To find an approximate solution, we discretize (2) using the Finite Element Method (FEM). Since we are operating in the $H(\text{curl})$ space, we construct our discretization using first-order Nédélec elements to generate basis functions $\varphi_i$. Letting $V_h$ be the space generated by our basis functions, we want to find $E_h \in V_h^g$ such that the variational formulation is satisfied for all $F_h \in V_h^0$, where $V_h^g$ denotes all elements of $V_h$ that adhere to boundary condition $g$. As in Example 1, the discretized variational formulation can be rewritten as a linear system $Ax=b$ because the basis functions are known. In particular, $A_{ij} = a(\varphi_i,\varphi_j)$, $b_i = b(\varphi_i)$, and $x$ is the vector of coefficients we are seeking in order to write our approximate solution as a linear combination of the basis functions.

In this example, first the mesh and finite elements are specified, then the linear form $b$ and bilinear form $A$ are assembled. The resulting linear system is then solved via preconditioned conjugate gradient iterations. Finally, the solution is recovered and the error is calculated against the manufactured solution $$E(x) = ( \sin(\kappa y), \sin(\kappa z), \sin(\kappa x)),$$ where $\kappa > 0$ is a prescribed frequency parameter.

  Annotated Example 3

The source file examples/ex3.cpp implements the electromagnetics problem described above. The main parts of the example code are annotated below.

Possible mesh options:

Many possible mesh options are written in lines 5-23 for the convenience of the user along with example runs for ex3.cpp with each of those mesh files. A selection of mesh files are available in MFEM directory in the ../mfem/data folder. The user can also specify a custom file.

Parse command-line options

We start the code (lines 70-94) with several command line options:

1. -m: Specifies the mesh file. MFEM can handle various mesh types. Please check the official documentation for further details.
2. -o: Specifies the polynomial degree of the finite element (Nédélec element in this case).
3. -f: Specifies the variable freq, which sets the frequency parameter kappa = freq * M_PI; in line 96; cf. (3).
4. -sc: Specifies whether or not to use Static Condensation in line 179.
5. -pa: Specifies whether or not to use Partial Assembly in line 171. This can reduce both the size of the problem and computational time.
6. -nc: Specifies whether or not the mesh should remain conforming during refinement.
7. -d: Specifies the computation device. Possible options are shown in lines 25-29. The user can choose to use RAJA, CUDA or, by default, cpu.
8. -vis: Determines if visualization will be executed or not (with the help of GLVis server in lines 229-236).

Device choice

   Device device(device_config);
   device.Print();

MFEM allows GPUs and other programming models such as CUDA, RAJA, OpenMP, and OCCA. The command line option -d provides the device to be used. In line 100, the device is set from the command-line input device_config. It is then printed in output to confirm the execution to the user in line line 101.

Mesh creation

   Mesh *mesh = new Mesh(mesh_file, 1, 1);
   dim = mesh->Dimension();
   int sdim = mesh->SpaceDimension();

This block of code defines the mesh after reading from the mesh_file provided by the user. The 1, 1 arguments control the generation of boundary elements and preparation for conforming or non-conforming refinement as selected by the -nc option. Then dim retrieves the dimension of the mesh elements (e.g., 2 for triangles, 3 for tetrahedrons). Then sdim retrieves the dimension of the ambient space.

Mesh refinement

   {
      int ref_levels =
         (int)floor(log(50000./mesh->GetNE())/log(2.)/dim);
      for (int l = 0; l < ref_levels; l++)
      {
         mesh->UniformRefinement();
      }
   }

This block of code refines the input mesh. It calculates the number of uniform refinement levels (ref_levels) needed so that the final mesh has roughly 50,000 elements or fewer. In line 121, the mesh refinement is done with a factor of 2 per dimension dim. Then in line 124, the mesh is refined ref_levels times in an uniform manner.

Finite element space definition

The next block of code defines the Finite Element Space on the mesh provided by the user.

   FiniteElementCollection *fec = new ND_FECollection(order, dim);
   FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
   cout << "Number of finite element unknowns: "
        << fespace->GetTrueVSize() << endl;

Here, the FiniteElementCollection corresponds to Nédélec basis functions because Nédélec finite elements have optimal approximation properties in $H(\mathrm{curl})$; i.e., the energy space over which problem (2) is well-defined.

Identification of essential boundary dofs

The next block of code marks the degrees of freedom associated with the conforming (i.e., 'true') essential/Dirichlet-type boundary conditions.

   Array<int> ess_tdof_list;
   if (mesh->bdr_attributes.Size())
   {
      Array<int> ess_bdr(mesh->bdr_attributes.Max());
      ess_bdr = 1;
      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
   }

Note that $E \times n$ is specified over $\partial\Omega$ in (2). Therefore, all boundary attributes are marked as essential by setting ess_bdr = 1 line 143.

Assembling the linear form

   VectorFunctionCoefficient f(sdim, f_exact);
   LinearForm *b = new LinearForm(fespace);
   b->AddDomainIntegrator(new VectorFEDomainLFIntegrator(f));
   b->Assemble();

This block of code defines the linear form on the right-hand side of (2). More specifically, VectorFEDomainLFIntegrator(f) encodes the action of the linear form $\psi_i \mapsto \int_{\partial\Omega} f \cdot \psi_i\, \mathrm{d}S$, where $\psi_i$ are the Nédélec basis functions. The corresponding right-hand side vector $b$ is assembled using the Assemble function.

Solution GridFunction

   GridFunction x(fespace);
   VectorFunctionCoefficient E(sdim, E_exact);
   x.ProjectCoefficient(E);

In this code block, the GridFunction x is constructed using fespace. This GridFunction will hold the dof values representing the finite element solution vector on the mesh. In order to match the Dirichlet boundary conditions, x is initialized by interpolating the exact solution E using the function ProjectCoefficient.

Defining the bilinear form

The bilinear form $a(E, F)$ in (2) is defined in this code block.

   Coefficient *muinv = new ConstantCoefficient(1.0);
   Coefficient *sigma = new ConstantCoefficient(1.0);
   BilinearForm *a = new BilinearForm(fespace);
   if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
   a->AddDomainIntegrator(new CurlCurlIntegrator(*muinv));
   a->AddDomainIntegrator(new VectorFEMassIntegrator(*sigma));

Here, muinv and sigma, which correspond to the two parameters appearing the strong formulation (1), are both set to 1.0. If partial assembly (pa) is enabled, it's set here in line 171. Then CurlCurlIntegrator is applied on a using the constant coefficient muinv, and subsequently, VectorFEMassIntegrator is applied on a using the constant coefficient sigma. Now we are ready to assemble the bilinear form in the below block.

Assembling the bilinear form

if (static_cond) { a->EnableStaticCondensation(); }
   a->Assemble();

   OperatorPtr A;
   Vector B, X;
   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);

   cout << "Size of linear system: " << A->Height() << endl;

Here, static condensation is enabled if it was requested from the command line using the command -sc. Then a->Assemble(); computes the element matrices and assembles the global stiffness matrix. The function FormLinearSystm transforms the assembled bilinear form a and linear form b into the final algebraic system A X = B; cf. lines 212-216 of Annotated Example 1). Note that the function FormLinearSystem also eliminates the essential boundary dofs declared in ess_tdof_list. Lastly, we print the size of the linear system for the convenience of the user.

Solving the linear system

The next block of code solves the linear system A X = B for the unknown dofs X. If Partial Assembly (pa) is used, it employs the Preconditioned Conjugate Gradient (PCG) iterative solver with a simple Jacobi preconditioner (OperatorJacobiSmoother). Otherwise, in case of full assembly, the code checks if the SuiteSparse library is used or not. If the SuiteSparse library has been used, then UMFPACK will be used to solve the system. Otherwise, a Gauss-Seidel preconditioner along with PCG will be used to solve the system.

if (pa) // Jacobi preconditioning in partial assembly mode
   {
      OperatorJacobiSmoother M(*a, ess_tdof_list);
      PCG(*A, M, B, X, 1, 1000, 1e-12, 0.0);
   }
   else
   {
#ifndef MFEM_USE_SUITESPARSE
      // 11. Define a simple symmetric Gauss-Seidel preconditioner and use it to
      //     solve the system Ax=b with PCG.
      GSSmoother M((SparseMatrix&)(*A));
      PCG(*A, M, B, X, 1, 500, 1e-12, 0.0);
#else
      // 11. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the
      //     system.
      UMFPackSolver umf_solver;
      umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
      umf_solver.SetOperator(*A);
      umf_solver.Mult(B, X);
#endif
   }

Observe that, PCG takes several arguments in line 200, which we now describe. A, B, X refers to the linear system A X = B. M refers to the solver after preconditioner applied (if applicable). 1 refers to the value of the variable print_iter, which implies a message about sending detailed information about each iteration. The maximum number of iterations allowed max_num_iter is 500. The relative tolerance allowed for convergence 1e-12 and the absolute tolerance allowed for convergence is 0.0. Users are encouraged to read the Solvers documentation for a better understanding.

Recovering the solution as a finite element grid function

After solving the system A X = B, the vector X only contains the values for the unknown (non-essential) degrees of freedom.

 a->RecoverFEMSolution(X, *b, x);

This line uses the RecoverFEMSolution method to reconstruct the complete solution GridFunction x. It copies the computed values from X into the corresponding entries of x, and adds back the known essential boundary DOF values.

Saving the refined mesh and the solution

This code block saves the final solution and refined mesh using line x.Save(sol_ofs);. It also prints the final mesh geometry. This mesh and solution can be visualized later using the GLVis server.

   {
    ofstream mesh_ofs("refined.mesh");
    mesh_ofs.precision(8);
    mesh->Print(mesh_ofs);
    ofstream sol_ofs("sol.gf");
    sol_ofs.precision(8);
    x.Save(sol_ofs);
   }

Visualization

If the visualization flag is enabled with the -vis command-line, the code block connects to a GLVis server and plots the results using the GLVis server. A new GLVis window opens which produces the customizable plot. For further details, please go through the GLVis Documentation.

   if (visualization)
   {
      char vishost[] = "localhost";
      int  visport   = 19916;
      socketstream sol_sock(vishost, visport);
      sol_sock.precision(8);
      sol_sock << "solution\n" << *mesh << x << flush;
   }

Freeing the used memory

Lastly, the variables are deleted to prevent a memory leak in C++, which might become an issue if the program is executed multiple times.

   delete a;
   delete sigma;
   delete muinv;
   delete b;
   delete fespace;
   delete fec;
   delete mesh;

Back to the  Getting Started page.