Many physical processes may be controlled by external constraints corresponding to complex boundary conditions in their theoretical models other than the standard Neumann, Dirichlet, or Robin boundary conditions. For example, in the electrically induced insulator-metal transitions Wook et al. (2008); Kim,Hyun-Tak et al. (2010); Shukla et al. (2014); Shi and Chen (2018, 2019); Kumar et al. (2017, 2020), the phase-changing material is often in series connection to a resistor, so a type of integral boundary condition naturally arises. The series resistor can protect the material from large current damage, and sometimes is even crucial for the emergence of desired phenomena such as the voltage oscillation in vanadium dioxide Wook et al. (2008); Kim,Hyun-Tak et al. (2010); Shukla et al. (2014) and the chaotic dynamics in niobium dioxide Mott memristors Kumar et al. (2017, 2020). This form of integral boundary condition may also arise in other physical processes. For example, one may encounter a scenario in which a material is in a finite thermal bath and they as a whole are isolated. These imply that the surface temperature of the material is uniform, and is determined by the energy conservation of the total energy of the material and thermal bath. The boundary condition on the surface of the material is then , where is the heat capacitance of the thermal bath, is the temperature on , is the time, and is the thermal conductivity of the material. In this paper, we discuss ways of handling such a special type of integral boundary condition encountered in many electric circuit systems. This boundary condition is illustrated by the typical circuit shown in Figure 1, in which a material of interest is electrically excited by a direct voltage source through a series resistor . We consider the material to be two-dimensional, that is, the material has uniform properties along the third dimension (out of plane). The thickness in the third dimension is set to . The length and width are both
. The physical partial differential equation (PDE) determining the electric potential fieldin the material is given by Gauss’s law with proper boundary conditions,
where denote the spatial coordinates, is the current passing through the boundary , is the fixed potential at the boundary , and the source term depends on the charge distribution inside the material. The special boundary condition arises at the boundary . The value of on is determined by Kirchhoff’s law and Ohm’s law, that is, it is the difference between the total voltage and the voltage drop across the resistor,
. We denote the area and the outward unit normal vector of the boundaryby and , respectively. The current is related to the current density (Ohm’s law, is the conductivity of the material) through the following surface integral:
The negative sign in the first equality arises since the current flows into the material across . Finally, we obtain the following integral boundary condition at :
The equipotential on is determined by an integral of its gradient over that boundary.
In formulating the variational problem for PDEs with such complex integral boundary conditions, and solving it by the finite element method, one difficulty is that these boundary conditions cannot be easily incorporated into the weak form like a Neumann or Robin boundary condition, or enforced on the boundaries like a Dirichlet boundary condition. Therefore, we have to develop special methods to deal with this kind of special boundary condition. A general way of handling a special boundary condition is to treat it as a constraint for the PDE, and then use the Lagrange multiplier method to convert the PDE problem to a stationary point problem. An example of treating a Dirichlet boundary condition using Lagrange multiplier method can be found in IBabuška (1973); Glowinski et al. (1995). Another widely used method is to add a penalty term such as the Nitsche’s method Nitsche (1971); Babuška (1973). In this paper, we propose a modified Nitsche’s method to handle the integral boundary condition (1.4) in solving the PDE (1.1) by the finite element method.
The paper is organized as follows. In Section 2, we formulate the variational problem using the Lagrange multiplier method, then we propose a modified Nitsche’s method that allows us to incorporate the special boundary condition into the variational form naturally. Numerical tests are presented in Section 3. The results show that our proposed method achieves better accuracy than the Lagrange multiplier method, and the resulting linear system can be efficiently solved by an iterative solver with algebraic multigrid (AMG) preconditioning. Finally we come to a conclusion in Section 4.
In this section, we briefly review the Lagrange multiplier method and then propose a new method to solve problem (1.1) with boundary condition (1.4), based on Nitsche’s method. In order to derive the variational form of the two methods, we define the following function spaces:
where , and .
For the space discretization, we use a uniform triangular mesh of the domain , where denotes the mesh size of . We let be the standard finite element subspaces of defined as follows,
where is the space of polynomials of degree at most on , with . In our numerical tests, we use continuous piecewise linear functions for the finite element space, i.e. .
2.1 Lagrange multiplier method
In this subsection, we use the Lagrange multiplier method to solve equation (1.1) with integral boundary condition (1.4) on . By this method ELSGOLC (1961), solving equation (1.1) with (1.4) is equivalent to finding the stationary point of the following functional
where , and is the Lagrange multiplier. Based on the functional (2.7), we introduce the following bilinear form
and a functional
where , and .
The stationary point , of (2.7) is such that for all , ,
The corresponding finite element formulation is shown in Algorithm 1.
2.2 Modified Nitsche’s method
In the Lagrange multiplier method (2.11), an additional variable is introduced. We also note that in the bilinear form of (2.11), there are derivative terms in the integrals, which might lead to convergence problems. Here we propose a modified Nitsche’s method to handle these issues. This proposed method does not require an additional Lagrange multiplier variable, and the resulting boundary condition does not involve a derivative term in the integral.
We assume that is always positive, i.e, , which is true in most physical situations. By introducing a small term with a small positive fine-tuned constant , we replace the boundary condition (1.4) by the following
Integrating both sides of equation (2.13) over and rearranging the equation, we obtain
which can be implemented like a Neumann boundary condition.
We introduce the following bilinear form
and a functional
Now the variational problem reads as follows: find , for all , such that
The well-posedness of the above variational problem is proved in Appendix for a special case when is a constant.
3 Numerical tests
In this section, we carry out numerical tests to compare the Lagrange multiplier method (Algorithm 1) with the modified Nitsche’s method (Algorithm 2). We implement the two methods using FEniCS Alnæs et al. (2015). Details of the weak form implementation can be found in Appendix. In the first test, is a constant, and in the second test, is a varying function. We first show that in both cases, our proposed modified Nitsche’s method in Algorithm 2 gives better accuracy. We further demonstrate that our proposed method can be solved using generalized minimal residual method (GMRES) Saad and Schultz (1986) with AMG preconditioning Stüben (1999).
Test 1. We consider a manufactured PDE in the same form as (1.1), with , , , , , . This manufactured PDE has a unique solution . The parameter used in the modified Nitsche’s method is set to . Numerical tests for a wide range of are shown in Appendix.
In this test, the resulting linear systems are solved by the default direct solver using Sparse LU decomposition. For a sequence of tests, the initial mesh size is set to and then refined by a factor of . For each mesh refinement, we report error inside the domain, , error on the boundary , , and error inside the domain, , i.e. , where is the exact solution and is the finite element solution. The results obtained by the Lagrange multiplier method and the modified Nitsche’s method are shown in Tables 1 and 2, respectively. Numerical solutions obtained by the modified Nitsche’s method achieve optimal convergence rate of nd order under the norm inside the domain. When the Lagrange multiplier method is used, the error convergence rate is not optimal, although the numerical solutions still converge. As to error on the boundary , the modified Nitsche’s method gives almost the exact solution whereas the numerical solution obtained by the Lagrange multiplier method is much less accurate. For error inside the domain, both methods achieve optimal convergence rate.
Test 2. We consider another manufactured PDE, with , , , , in . This manufactured PDE has a unique solution . The parameter is also set to . The testing results for the Lagrange multiplier method and the modified Nitsche’s method are shown in Tables 3 and 4, respectively. The numerical solutions obtained by the modified Nitsche’s method again achieve optimal convergence rate of nd order under the norm inside the domain. For solutions obtained by the Lagrange multiplier method, the error convergence rate is still not optimal. As to error on the boundary , the modified Nitsche’s method also achieves nd order convergence rate while the Lagrange multiplier method only has 1st order convergence rate.
We further investigate the linear solver GMRES with AMG preconditioning for solving the manufactured PDEs in both Test 1 and Test 2. We record in Tables 5 and 6 the number of iterations needed for convergence in Tests 1 and 2, respectively. The stopping criterion is when the actual residual norm is smaller than . It can be seen that the iterative solver GMRES with AMG preconditioning can be directly applied to the linear system resulted from the modified Nitsche’s method. Moreover, the number of iterations does not increase much when the mesh is refined, which indicates the effectiveness of AMG preconditioning for the modified Nitsche’s method. However, under the same settings, the Lagrange multiplier method does not converge. The ability of being solved with AMG preconditioning is also an advantage of our modified Nitsche’s method over the Lagrange multiplier method, which can be crucially benefiting in addressing large scale systems.
|Mesh size||Modified Nitsche’s method||Lagrange multiplier method|
|Mesh size||Modified Nitsche’s method||Lagrange multiplier method|
We derived a special integral boundary condition for the Poisson equation for a typical electric circuit model. We proposed a modified Nitsche’s method to deal with the integral boundary condition and compared it with the Lagrange multiplier method. The proposed method computed the numerical solutions more accurately than the Lagrange multiplier method did, and achieved optimal convergence rate under both and norms. Moreover, iterative solvers equipped with AMG preconditioner can be direcly used for solving the linear system resulted from our modified Nitsche’s method, which however did not converge for the Lagrange multiplier method. The modified Nitsche’s method can be useful in many physical models possessing this kind of integral boundary condition.
This work is supported as part of the Computational Materials Sciences Program funded by the US Department of Energy, Office of Science, Basic Energy Sciences, under Award Number DE-SC0020145.
Appendix A. Implementation of the Lagrange multiplier method in FEniCS
The resulting variational form (2.10) can not be directly implemented using FEniCS due to the two double integrals in (2.1). Hence we introduce an auxiliary variable and reformulate the problem. It then becomes seeking for the stationary point of the following functional
under the constraint that
Let be the space of real numbers, and let be the test function of . The variational form consists of two parts. From (A.1), we have
From (A.2), we can derive the following weak form,
where is the length of the boundary .
Appendix B. Implementation of the modified Nitsche’s method in FEniCS
Due to the same reason as in the Lagrange multiplier method, in the actual implementation, we introduce an auxiliary variable for the variational form (2.19). Again we can derive the variational form for the auxiliary variable, which is
Appendix C. Test of
We look at how different choices of affect the numerical accuracy. For this test, we use the manufactured PDE in Test 2 and compute the numerical solution on the mesh with size .
We can see the modified Nitsche’s method works for a wide range of given a certain error tolerance. Also, as expected, the accuracy is poor with big . This is because the modified boundary condition (2.16) becomes a poor approximation to the original integral boundary condition (1.4) when is large. In the case when is very small, from to the numerical accuracy also deteriorates.
Appendix D. Well-posedness
We prove the well-posedness of the variational problem (2.19) for a special case when is a constant.
Assuming that ,the variational problem (2.19) for the modified Nitsche’s method is well-posed for any if is a constant.
Without loss of generality, we may assume . So and are the same function space. Since is a constant, we can write the bilinear form in (2.19) as follows,
To see its boundedness property, we have
where the last inequality follows from the trace theorem.
Knowing is a constant and applying Cauchy Schwarz inequality , where is the length of , we have
So is positive definite. By Poincare inequality, we have that the bilinear form is strictly positive definite. Since is strictly positive definite and is positive definite, the bilinear form must be strictly positive definite.
Furthermore, following a similar argument for the bounded property of , we have
This shows that is a bounded linear functional on .
Therefore, the well-posedness follows from the Lax-Milgram lemma Bressan (2013). ∎
6 Data Availability
- The fenics project version 1.5. Archive of Numerical Software 3 (100). External Links: Cited by: §3.
- The finite element method with penalty. Mathematics of computation 27 (122), pp. 221–228. Cited by: §1.
- Lecture notes on functional analysis: with applications to linear partial differential equations. Graduate studies in mathematics, American Mathematical Society. External Links: Cited by: Appendix D. Well-posedness.
- CHAPTER i - the method of variation in problems with fixed boundaries. In Calculus of Variations, L.E. ELSGOLC (Ed.), International Series of Monographs on Pure and Applied Mathematics, pp. 13 – 63. External Links: Cited by: §2.1.
- A lagrange multiplier/fictitious domain method for the dirichlet problem—generalization to some flow problems. Japan journal of industrial and applied mathematics 12 (1), pp. 87. Cited by: §1.
- The finite element method with lagrangian multipliers. Numerische Mathematik 20 (3), pp. 179–192. External Links: Cited by: §1.
- Electrical oscillations induced by the metal-insulator transition in vo2. Journal of Applied Physics 107 (2), pp. 023702. External Links: Cited by: §1.
- Chaotic dynamics in nanoscale mott memristors for analogue computing. Nature 548 (7667), pp. 318–321. External Links: Cited by: §1.
- Third-order nanocircuit elements for neuromorphic engineering. Nature 585 (7826), pp. 518–523. External Links: Cited by: §1.
- Über ein variationsprinzip zur lösung von dirichlet-problemen bei verwendung von teilräumen, die keinen randbedingungen unterworfen sind. Abhandlungen aus dem Mathematischen Seminar der Universität Hamburg 36 (1), pp. 9–15. External Links: Cited by: §1, §2.2.
- GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM Journal on Scientific and Statistical Computing 7 (3), pp. 856–869. External Links: Cited by: §3.
- Phase-field model of insulator-to-metal transition in under an electric field. Phys. Rev. Materials 2 (5), pp. 053803. External Links: Cited by: §1.
- Current-driven insulator-to-metal transition in strongly correlated . Phys. Rev. Applied 11, pp. 014059. External Links: Cited by: §1.
- Synchronized charge oscillations in correlated electron systems. Scientific reports 4, pp. 4964. External Links: Cited by: §1.
- Algebraic Multigrid (AMG): An Introduction with Applications. GMD-Forschungszentrum Informationstechnik. Cited by: §3.
- Metal-insulator transition-induced electrical oscillation in vanadium dioxide thin film. Applied Physics Letters 92 (16), pp. 162903. External Links: Cited by: §1.