# Integral boundary conditions in phase field models

Xiaofeng Xu  Lian Zhang  Yin Shi  Long-Qing Chen  Jinchao Xu 111Corresponding
author<EMAIL_ADDRESS>(Jinchao Xu) Department of Mathematics, Pennsylvania
State University, University Park, PA, 16802, USA  Department of Material
Science and Engineering, Pennsylvania State University, University Park, PA,
16802, USA

###### Abstract

Modeling the microstructure evolution of a material embedded in a device often
involves integral boundary conditions. Here we propose a modified Nitsche’s
method to solve the Poisson equation with an integral boundary condition,
which is coupled to phase-field equations of the microstructure evolution of a
strongly correlated material undergoing metal-insulator transitions. Our
numerical experiments demonstrate that the proposed method achieves optimal
convergence rate while the rate of convergence of the conventional Lagrange
multiplier method is not optimal. Furthermore, the linear system derived from
the modified Nitsche’s method can be solved by an iterative solver with
algebraic multigrid preconditioning. The modified Nitsche’s method can be
applied to other physical boundary conditions mathematically similar to this
electric integral boundary condition.

###### keywords:

Phase field, Integral boundary condition, Lagrange multiplier method, Modified
Nitsche’s method, Algebraic multigrid.

##  1 Introduction

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 [1, 2, 3, 4,
5, 6, 7], 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 [1, 2, 3] and the chaotic dynamics in
niobium dioxide Mott memristors [6, 7]. 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
$S$ of the material is then $C_{b}\partial T/\partial t=-\int_{S}k(\partial
T/\partial n)dA$, where $C_{b}$ is the heat capacitance of the thermal bath,
$T$ is the temperature on $S$, $t$ is the time, and $k$ 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 $V$ through a series resistor $R$. 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 $1$. The length and width are both $L$. The physical partial
differential equation (PDE) determining the electric potential field
$\Phi(x,y)$ in the material is given by Gauss’s law with proper boundary
conditions,

$\left\\{\begin{aligned} -\Delta\Phi(x,y)&=f(x,y)&&\hbox{~{}~{}~{}~{}in
$\Omega=(0,L)^{2}$},\\\ \Phi&=V-IR&&\hbox{~{}~{}~{}~{}on
$\Gamma_{1}=\\{(x,y)\in\partial\Omega:x=0\\}$},\\\
\frac{\partial\Phi}{\partial y}&=0&&\hbox{~{}~{}~{}~{}on
$\Gamma_{2}=\\{(x,y)\in\partial\Omega:y=0~{}\text{or}~{}L\\}$},\\\
\Phi&=\Phi_{D}&&\hbox{~{}~{}~{}~{}on
$\Gamma_{3}=\\{(x,y)\in\partial\Omega:x=L\\}$},\end{aligned}\right.$ (1.1)

where $(x,y)$ denote the spatial coordinates, $I$ is the current passing
through the boundary $\Gamma_{1}$, $\Phi_{D}$ is the fixed potential at the
boundary $\Gamma_{3}$, and the source term $f(x,y)$ depends on the charge
distribution inside the material. The special boundary condition arises at the
boundary $\Gamma_{1}$. The value of $\Phi$ on $\Gamma_{1}$ is determined by
Kirchhoff’s law and Ohm’s law, that is, it is the difference between the total
voltage $V$ and the voltage drop across the resistor, $IR$. We denote the area
and the outward unit normal vector of the boundary $\Gamma_{1}$ by
$S=\Gamma_{1}\times[0,1]$ and $\bm{n}$, respectively. The current $I$ is
related to the current density $\bm{J}=-\sigma\nabla\Phi$ (Ohm’s law, $\sigma$
is the conductivity of the material) through the following surface integral:

$\displaystyle I$ $\displaystyle=-\int_{S}\bm{J}\cdot\bm{n}dA$ (1.2)
$\displaystyle=\int_{\Gamma_{1}}\sigma\frac{\partial\Phi}{\partial n}ds.$
(1.3)

The negative sign in the first equality arises since the current flows into
the material across $\Gamma_{1}$. Finally, we obtain the following integral
boundary condition at $\Gamma_{1}$:

$\displaystyle\Phi=V-R\int_{\Gamma_{1}}\sigma\frac{\partial\Phi}{\partial
n}ds.$ (1.4)

The equipotential on $\Gamma_{1}$ is determined by an integral of its gradient
over that boundary.

Figure 1: A simple model

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
[8, 9]. Another widely used method is to add a penalty term such as the
Nitsche’s method [10, 11]. 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.

##  2 Methods

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:

$\displaystyle V_{\Phi}$ $\displaystyle=\left\\{\Phi\in
H^{1}({\Omega})~{}|~{}\Phi|_{\Gamma_{3}}=\Phi_{D}\right\\},$ (2.1)
$\displaystyle V_{\Phi,0}$ $\displaystyle=\left\\{\Phi\in
H^{1}({\Omega})~{}|~{}\Phi|_{\Gamma_{3}}=0\right\\},$ (2.2) $\displaystyle
V_{\lambda}$ $\displaystyle=L^{2}(\Gamma_{1}),$ (2.3)

where
$L^{2}(\Gamma_{1})=\left\\{\Phi:\Gamma_{1}\to\mathbb{R}~{}|~{}\int_{\Gamma_{1}}\Phi^{2}dx<\infty\right\\}$,
and $H^{1}({\Omega})=\left\\{\Phi\in
L^{2}(\Omega)~{}|~{}\int_{\Omega}|\nabla\Phi|^{2}dx<\infty\right\\}$.

For the space discretization, we use a uniform triangular mesh
$\mathcal{T}^{h}$ of the domain $\Omega$, where $h$ denotes the mesh size of
$\mathcal{T}^{h}$. We let $V^{h}_{\Phi},V^{h}_{\Phi,0},V^{h}_{\lambda}$ be the
standard finite element subspaces of $V_{\Phi},V_{\Phi,0},V_{\lambda}$ defined
as follows,

$\displaystyle V^{h}_{\Phi}$ $\displaystyle=\left\\{\Phi_{h}\in
C(\overline{\Omega})\cap V_{\Phi}~{}|~{}\Phi_{h}|_{\mathcal{K}}\in
P_{n}(\mathcal{K}),\forall\mathcal{K}\in\mathcal{T}^{h}\right\\},$ (2.4)
$\displaystyle V^{h}_{\Phi,0}$ $\displaystyle=\left\\{\Phi_{h}\in
C(\overline{\Omega})\cap V_{\Phi,0}~{}|~{}\Phi_{h}|_{\mathcal{K}}\in
P_{n}(\mathcal{K}),\forall\mathcal{K}\in\mathcal{T}^{h}\right\\},$ (2.5)
$\displaystyle V^{h}_{\lambda}$ $\displaystyle=\left\\{\lambda_{h}\in
C(\Gamma_{1})~{}|~{}\lambda_{h}|_{\mathcal{K}\cap\Gamma_{1}}\in
P_{n}(\mathcal{K}\cap\Gamma_{1}),\forall\mathcal{K}\in\mathcal{T}^{h}\right\\},$
(2.6)

where $P_{n}(\mathcal{K})$ is the space of polynomials of degree at most $n$
on $\mathcal{K}\in\mathcal{T}^{h}$, with $n\geq 1$. In our numerical tests, we
use continuous piecewise linear functions for the finite element space, i.e.
$n=1$.

###  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 $\Gamma_{1}$. By this method
[12], solving equation (1.1) with (1.4) is equivalent to finding the
stationary point of the following functional

$\displaystyle\int_{\Omega}\left(\frac{1}{2}|\nabla\Phi|^{2}-f\Phi\right)dx-\int_{\Gamma_{1}}\lambda\left(V-\Phi-R\int_{\Gamma_{1}}\sigma\frac{\partial\Phi}{\partial
n}ds\right)ds,$ (2.7)

where $\Phi\in V_{\Phi}$, and $\lambda\in V_{\lambda}$ is the Lagrange
multiplier. Based on the functional (2.7), we introduce the following bilinear
form

$\displaystyle B(\Phi,\lambda;v_{1},v_{2})$
$\displaystyle=\int_{\Omega}\nabla\Phi\cdot\nabla
v_{1}dx+\int_{\Gamma_{1}}\lambda v_{1}ds$
$\displaystyle+R\int_{\Gamma_{1}}\lambda\bigg{(}\int_{\Gamma_{1}}\sigma\frac{\partial
v_{1}}{\partial n}ds\bigg{)}ds+\int_{\Gamma_{1}}\Phi
v_{2}ds+R\int_{\Gamma_{1}}\bigg{(}\int_{\Gamma_{1}}\sigma\frac{\partial\Phi}{\partial
n}ds\bigg{)}v_{2}ds,$ (2.8)

and a functional

$\displaystyle
F(v_{1},v_{2})=\int_{\Omega}fv_{1}dx+V\int_{\Gamma_{1}}v_{2}ds,$ (2.9)

where $v_{1}\in V_{\Phi,0}$, and $v_{2}\in V_{\lambda}$.

The stationary point $(\Phi$, $\lambda)\in V_{\Phi}\times V_{\lambda}$ of
(2.7) is such that for all $v_{1}\in V_{\Phi,0}$, $v_{2}\in V_{\lambda}$,

$\displaystyle B(\Phi,\lambda;v_{1},v_{2})=F(v_{1},v_{2}).$ (2.10)

The corresponding finite element formulation is shown in Algorithm 1.

Find $\Phi_{h}\in V_{\Phi}^{h}$, $\lambda_{h}\in V_{\lambda}^{h}$ such that
for all $v_{1,h}\in V_{\Phi,0}^{h}$, $v_{2,h}\in V_{\lambda}^{h}$,
$\displaystyle B(\Phi_{h},\lambda_{h};v_{1,h},v_{2,h})=F(v_{1,h},v_{2,h}).$
(2.11)

Algorithm 1 FEM with Lagrange multiplier method

###  2.2 Modified Nitsche’s method

In the Lagrange multiplier method (2.11), an additional variable $\lambda_{h}$
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 $\sigma$ is always positive, i.e,
$\sigma:[0,L]\to\mathbb{R}^{+}$, which is true in most physical situations. By
introducing a small term $\epsilon\frac{\partial\Phi}{\partial n}$ with a
small positive fine-tuned constant $\epsilon$, we replace the boundary
condition (1.4) by the following

$\displaystyle\Phi+\epsilon\frac{\partial\Phi}{\partial
n}=V-R\int_{\Gamma_{1}}\sigma\frac{\partial\Phi}{\partial n}ds.$ (2.12)

Using (2.12) directly as a boundary condition is the classic Nitsche’s method
[10]. Multiplying $\sigma$ on both sides of equation (2.12), we have

$\displaystyle\sigma\Phi+\epsilon\sigma\frac{\partial\Phi}{\partial n}=\sigma
V-R\sigma\int_{\Gamma_{1}}\sigma\frac{\partial\Phi}{\partial n}ds.$ (2.13)

Integrating both sides of equation (2.13) over $\Gamma_{1}$ and rearranging
the equation, we obtain

$\displaystyle\int_{\Gamma_{1}}\sigma\frac{\partial\Phi}{\partial
n}ds=\frac{1}{R\int_{\Gamma_{1}}\sigma ds+\epsilon}\int_{\Gamma_{1}}(\sigma
V-\sigma\Phi)ds.$ (2.14)

Substituting equation (2.14) into equation (2.12) gives

$\displaystyle\Phi+\epsilon\frac{\partial\Phi}{\partial
n}=V-\frac{R}{R\int_{\Gamma_{1}}\sigma ds+\epsilon}\int_{\Gamma_{1}}(\sigma
V-\sigma\Phi)dy.$ (2.15)

Rearranging equation (2.15), we obtain the boundary condition on $\Gamma_{1}$
for problem (1.1)

$\displaystyle\frac{\partial\Phi}{\partial
n}=\frac{1}{\epsilon}\bigg{(}V-\Phi-\frac{R}{R\int_{\Gamma_{1}}\sigma
ds+\epsilon}\int_{\Gamma_{1}}\sigma Vds+\frac{R}{R\int_{\Gamma_{1}}\sigma
ds+\epsilon}\int_{\Gamma_{1}}\sigma\Phi ds\bigg{)},$ (2.16)

which can be implemented like a Neumann boundary condition.

We introduce the following bilinear form

$a(\Phi,v_{1})=\int_{\Omega}\nabla\Phi\cdot\nabla
v_{1}dx+\frac{1}{\epsilon}\int_{\Gamma_{1}}\bigg{(}\Phi-\frac{R}{R\int_{\Gamma_{1}}\sigma
ds+\epsilon}\int_{\Gamma_{1}}\sigma\Phi ds\bigg{)}v_{1}ds,$ (2.17)

and a functional

$L(v_{1})=\int_{\Omega}fv_{1}dx+\frac{1}{\epsilon}\int_{\Gamma_{1}}\bigg{(}V-\frac{R}{R\int_{\Gamma_{1}}\sigma
ds+\epsilon}\int_{\Gamma_{1}}\sigma Vds\bigg{)}v_{1}ds.$ (2.18)

Now the variational problem reads as follows: find $\Phi\in V_{\Phi}$, for all
$v_{1}\in V_{\Phi,0}$, such that

$a(\Phi,v_{1})=L(v_{1}).$ (2.19)

The well-posedness of the above variational problem is proved in Appendix for
a special case when $\sigma$ is a constant.

We discretize the variational problem (2.19) using FEM, which gives Algorithm
2.

Find $\Phi_{h}\in V_{\Phi}^{h}$ such that for all $v_{1,h}\in V_{\Phi,0}^{h}$,
$\displaystyle a(\Phi_{h},v_{1,h})=L(v_{1,h})$ (2.20)

Algorithm 2 FEM with modified Nitsche’s method.

##  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 [13]. Details of the weak form
implementation can be found in Appendix. In the first test, $\sigma$ is a
constant, and in the second test, $\sigma$ 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) [14] with AMG
preconditioning [15].

Test 1. We consider a manufactured PDE in the same form as (1.1), with
$f(x,y)=-4xy+2x$, $\sigma=1$, $V=1$, $R=1$,
$\Phi_{D}=\frac{2}{3}y^{3}-y^{2}+\frac{5}{6}$, $\Omega=(0,1)^{2}$. This
manufactured PDE has a unique solution
$\Phi_{e}=\frac{2}{3}xy^{3}-xy^{2}+\frac{5}{6}$. The parameter $\epsilon$ used
in the modified Nitsche’s method is set to $10^{-9}$. Numerical tests for a
wide range of $\epsilon$ 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 $0.1$ and then refined by a factor of $2$. For each mesh
refinement, we report $L^{2}$ error inside the domain,
$\|\Phi_{e}-\Phi_{h}\|_{L^{2}(\Omega)}$, $L^{2}$ error on the boundary
$\Gamma_{1}$, $\|\Phi_{e}-\Phi_{h}\|_{L^{2}(\Gamma_{1})}$, and $H^{1}_{0}$
error inside the domain, $\|\Phi_{e}-\Phi_{h}\|_{H^{1}_{0}(\Omega)}$, i.e.
$\|\nabla(\Phi_{e}-\Phi_{h})\|_{L^{2}(\Omega)}$, where $\Phi_{e}$ is the exact
solution and $\Phi_{h}$ 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 $2$nd order under the
$L^{2}$ norm inside the domain. When the Lagrange multiplier method is used,
the $L^{2}$ error convergence rate is not optimal, although the numerical
solutions still converge. As to $L^{2}$ error on the boundary $\Gamma_{1}$,
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 $H^{1}_{0}$ error inside the domain, both methods achieve
optimal convergence rate.

Mesh size | $\|\Phi_{e}-\Phi_{h}\|_{L^{2}(\Omega)}$ | order | $\|\Phi_{e}-\Phi_{h}\|_{H^{1}_{0}(\Omega)}$ | order | $\|\Phi_{e}-\Phi_{h}\|_{L^{2}(\Gamma_{1})}$ | order
---|---|---|---|---|---|---
$0.1$ | $5.10\times 10^{-3}$ | * | $2.96\times 10^{-2}$ | * | $8.33\times 10^{-3}$ | *
$0.05$ | $2.47\times 10^{-3}$ | 1.04 | $1.49\times 10^{-2}$ | 1.00 | $4.17\times 10^{-3}$ | 1.00
$0.025$ | $1.22\times 10^{-3}$ | 1.02 | $7.43\times 10^{-3}$ | 1.00 | $2.08\times 10^{-3}$ | 1.00
$0.0125$ | $6.05\times 10^{-4}$ | 1.01 | $3.72\times 10^{-3}$ | 1.00 | $1.04\times 10^{-3}$ | 1.00

Table 1:  Test 1 Error computation: Lagrange multiplier method with the
default direct solver

Mesh size | $\|\Phi_{e}-\Phi_{h}\|_{L^{2}(\Omega)}$ | order | $\|\Phi_{e}-\Phi_{h}\|_{H^{1}_{0}(\Omega)}$ | order | $\|\Phi_{e}-\Phi_{h}\|_{L^{2}(\Gamma_{1})}$ | order
---|---|---|---|---|---|---
$0.1$ | $7.35\times 10^{-4}$ | * | $2.84\times 10^{-2}$ | * | $2.01\times 10^{-7}$ | *
$0.05$ | $1.85\times 10^{-4}$ | 1.99 | $1.43\times 10^{-2}$ | 1.00 | $9.88\times 10^{-8}$ | *
$0.025$ | $4.63\times 10^{-5}$ | 2.00 | $7.13\times 10^{-3}$ | 1.00 | $1.89\times 10^{-7}$ | *
$0.0125$ | $1.16\times 10^{-5}$ | 2.00 | $3.57\times 10^{-3}$ | 1.00 | $1.13\times 10^{-7}$ | *

Table 2: Test 1 Error computation: modified Nitsche’s method with the default
direct solver

Test 2. We consider another manufactured PDE, with $f(x,y)=\sin x\cos\pi
y+\pi^{2}\sin x\cos\pi y$, $\sigma=y+1$, $V=1+\frac{2}{\pi^{2}}$, $R=1$,
$\Phi_{D}=1+\sin(1)\cos\pi y$ in $\Omega=(0,1)^{2}$. This manufactured PDE has
a unique solution $\Phi_{e}=\sin x\cos\pi y+1$. The parameter $\epsilon$ is
also set to $10^{-9}$. 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 $2$nd order under the $L^{2}$ norm inside
the domain. For solutions obtained by the Lagrange multiplier method, the
$L^{2}$ error convergence rate is still not optimal. As to $L^{2}$ error on
the boundary $\Gamma_{1}$, the modified Nitsche’s method also achieves $2$nd
order convergence rate while the Lagrange multiplier method only has 1st order
convergence rate.

Mesh size | $\|\Phi_{e}-\Phi_{h}\|_{L^{2}(\Omega)}$ | order | $\|\Phi_{e}-\Phi_{h}\|_{H^{1}_{0}(\Omega)}$ | order | $\|\Phi_{e}-\Phi_{h}\|_{L^{2}(\Gamma_{1})}$ | order
---|---|---|---|---|---|---
$0.1$ | $3.51\times 10^{-2}$ | * | $1.62\times 10^{-1}$ | * | $5.81\times 10^{-2}$ | *
$0.05$ | $1.74\times 10^{-2}$ | 1.01 | $8.16\times 10^{-2}$ | 0.99 | $2.96\times 10^{-2}$ | 0.97
$0.025$ | $8.69\times 10^{-3}$ | 1.00 | $4.08\times 10^{-2}$ | 1.00 | $1.49\times 10^{-2}$ | 0.99
$0.0125$ | $4.34\times 10^{-3}$ | 1.00 | $2.04\times 10^{-2}$ | 1.00 | $7.47\times 10^{-3}$ | 0.99

Table 3:  Test 2 Error computation: Lagrange multiplier method with the
default direct solver

Mesh size | $\|\Phi_{e}-\Phi_{h}\|_{L^{2}(\Omega)}$ | order | $\|\Phi_{e}-\Phi_{h}\|_{H^{1}_{0}(\Omega)}$ | order | $\|\Phi_{e}-\Phi_{h}\|_{L^{2}(\Gamma_{1})}$ | order
---|---|---|---|---|---|---
$0.1$ | $4.85\times 10^{-3}$ | * | $1.52\times 10^{-1}$ | * | $1.91\times 10^{-3}$ | *
$0.05$ | $1.22\times 10^{-3}$ | 1.99 | $7.60\times 10^{-2}$ | 1.00 | $4.77\times 10^{-4}$ | 2.00
$0.025$ | $3.05\times 10^{-4}$ | 2.00 | $3.80\times 10^{-2}$ | 1.00 | $1.19\times 10^{-4}$ | 2.00
$0.0125$ | $7.63\times 10^{-5}$ | 2.00 | $1.90\times 10^{-2}$ | 1.00 | $3.02\times 10^{-5}$ | 1.98

Table 4: Test 2 Error computation: modified Nitsche’s method with the default
direct solver

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 $10^{-7}$. 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
---|---|---
$0.1$ | 7 | Not Converge
$0.05$ | 7 | Not Converge
$0.025$ | 7 | Not Converge
$0.0125$ | 8 | Not Converge

Table 5: Number of iterations needed using GMRES and AMG preconditioning for
the PDE in Test 1

Mesh size | Modified Nitsche’s method | Lagrange multiplier method
---|---|---
$0.1$ | 7 | Not Converge
$0.05$ | 8 | Not Converge
$0.025$ | 9 | Not Converge
$0.0125$ | 10 | Not Converge

Table 6: Number of iterations needed using GMRES and AMG preconditioning for
the PDE in Test 2

##  4 Conclusion

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 $L^{2}$ and $H^{1}_{0}$ 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.

##  5 Acknowledgement

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 $w=\int_{\Gamma_{1}}\sigma\frac{\partial\Phi}{\partial
n}ds$ and reformulate the problem. It then becomes seeking for the stationary
point of the following functional

$\displaystyle\int_{\Omega}\big{(}\frac{1}{2}|\nabla\Phi|^{2}-f\Phi\big{)}dx-\int_{\Gamma_{1}}\lambda\big{(}V-\Phi-
Rw\big{)}ds,$ (A.1)

under the constraint that

$\displaystyle w=\int_{\Gamma_{1}}\sigma\frac{\partial\Phi}{\partial n}ds.$
(A.2)

Let $V_{3}$ be the space of real numbers, and let $v_{3}\in V_{3}$ be the test
function of $w$. The variational form consists of two parts. From (A.1), we
have

$\displaystyle\int_{\Omega}(\nabla\Phi\cdot\nabla
v_{1}-fv_{1})dx+\int_{\Gamma_{1}}\lambda
v_{1}ds-\int_{\Gamma_{1}}v_{2}(V-\Phi-Rw)ds+R\int_{\Gamma_{1}}\lambda
v_{3}ds=0$ (A.3)

From (A.2), we can derive the following weak form,

$\displaystyle\int_{\Gamma_{1}}\big{(}w-\int_{\Gamma_{1}}\sigma\frac{\partial\Phi}{\partial
n}ds\big{)}v_{3}ds=0$ (A.4) $\displaystyle\iff$
$\displaystyle\int_{\Gamma_{1}}wv_{3}ds-\int_{\Gamma_{1}}\sigma\frac{\partial\Phi}{\partial
n}ds\int_{\Gamma_{1}}v_{3}ds=0$ (A.5) $\displaystyle\iff$
$\displaystyle\int_{\Gamma_{1}}\big{(}w-|\Gamma_{1}|\sigma\frac{\partial\Phi}{\partial
n}\big{)}v_{3}ds=0,$ (A.6)

where $|\Gamma_{1}|$ is the length of the boundary $\Gamma_{1}$.

Combining (A.3) and (A.6), we obtain the variational form in the actual
implementation,

$\displaystyle\int_{\Omega}\nabla\Phi\cdot\nabla
v_{1}dx+\int_{\Gamma_{1}}\lambda
v_{1}ds+\int_{\Gamma_{1}}(\Phi+Rw)v_{2}ds+R\int_{\Gamma_{1}}\lambda
v_{3}ds+\int_{\Gamma_{1}}\big{(}w-|\Gamma_{1}|\sigma\frac{\partial\Phi}{\partial
n}\big{)}v_{3}ds=\int_{\Omega}fv_{1}dx+\int_{\Gamma_{1}}Vv_{2}ds$ (A.7)

## 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
$w=\int_{\Gamma_{1}}\sigma\Phi ds$ for the variational form (2.19). Again we
can derive the variational form for the auxiliary variable, which is

$\displaystyle\int_{\Gamma_{1}}\big{(}w-|\Gamma_{1}|\sigma\Phi\big{)}v_{3}ds=0,$
(B.1)

Combining (2.19) and (B.1), we obtain the following variational form in the
actual implementation,

$\displaystyle\int_{\Omega}\nabla\Phi\cdot\nabla
v_{1}dx+\frac{1}{\epsilon}\int_{\Gamma_{1}}\bigg{(}\Phi-\frac{R}{R\int_{\Gamma_{1}}\sigma
ds+\epsilon}w\bigg{)}v_{1}ds+\int_{\Gamma_{1}}\big{(}w-|\Gamma_{1}|\sigma\Phi\big{)}v_{3}ds=$
(B.2)
$\displaystyle\int_{\Omega}fv_{1}dx+\frac{1}{\epsilon}\int_{\Gamma_{1}}\bigg{(}V-\frac{R}{R\int_{\Gamma_{1}}\sigma
ds+\epsilon}\int_{\Gamma_{1}}\sigma Vds\bigg{)}v_{1}ds$ (B.3)

## Appendix C. Test of $\epsilon$

We look at how different choices of $\epsilon$ 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 $0.0125$.

$\epsilon$ | $\|\Phi_{e}-\Phi_{h}\|_{L^{2}(\Omega)}$ | $\|\Phi_{e}-\Phi_{h}\|_{H^{1}_{0}(\Omega)}$ | $\|\Phi_{e}-\Phi_{h}\|_{L^{2}(\Gamma_{1})}$
---|---|---|---
$10^{-2}$ | $3.03\times 10^{-3}$ | $2.26\times 10^{-2}$ | $7.28\times 10^{-3}$
$10^{-3}$ | $2.88\times 10^{-4}$ | $1.90\times 10^{-2}$ | $7.40\times 10^{-4}$
$10^{-4}$ | $6.76\times 10^{-5}$ | $1.90\times 10^{-2}$ | $7.08\times 10^{-5}$
$10^{-5}$ | $7.48\times 10^{-5}$ | $1.90\times 10^{-2}$ | $2.81\times 10^{-5}$
$10^{-6}$ | $7.60\times 10^{-5}$ | $1.90\times 10^{-2}$ | $2.95\times 10^{-5}$
$10^{-7}$ | $7.62\times 10^{-5}$ | $1.90\times 10^{-2}$ | $2.98\times 10^{-5}$
$10^{-8}$ | $7.62\times 10^{-5}$ | $1.90\times 10^{-2}$ | $2.98\times 10^{-5}$
$10^{-9}$ | $7.63\times 10^{-5}$ | $1.90\times 10^{-2}$ | $3.02\times 10^{-5}$
$10^{-10}$ | $7.70\times 10^{-5}$ | $1.90\times 10^{-2}$ | $3.26\times 10^{-5}$
$10^{-11}$ | $9.66\times 10^{-5}$ | $1.90\times 10^{-2}$ | $8.61\times 10^{-5}$
$10^{-12}$ | $2.85\times 10^{-4}$ | $1.90\times 10^{-2}$ | $4.45\times 10^{-4}$
$10^{-13}$ | $3.52\times 10^{-3}$ | $2.00\times 10^{-2}$ | $6.07\times 10^{-3}$
$10^{-14}$ | $2.72\times 10^{-2}$ | $5.08\times 10^{-2}$ | $4.71\times 10^{-2}$

Table C.1: Error computation with different epsilon values

We can see the modified Nitsche’s method works for a wide range of $\epsilon$
given a certain error tolerance. Also, as expected, the accuracy is poor with
big $\epsilon$. This is because the modified boundary condition (2.16) becomes
a poor approximation to the original integral boundary condition (1.4) when
$\epsilon$ is large. In the case when $\epsilon$ is very small, from
$10^{-14}$ to $10^{-11},$ 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 $\sigma$ is a constant.

######  Theorem 1.

Assuming that $f\in L^{2}(\Omega)$,the variational problem (2.19) for the
modified Nitsche’s method is well-posed for any $\epsilon>0$ if $\sigma$ is a
constant.

###### Proof.

Without loss of generality, we may assume $\Phi_{D}=0$. So $V_{\Phi}$ and
$V_{\Phi,0}$ are the same function space. Since $\sigma$ is a constant, we can
write the bilinear form $a(\cdot,\cdot)$ in (2.19) as follows,

$\displaystyle a(\Phi,v_{1})=a_{1}(\Phi,v_{1})+a_{2}(\Phi,v_{1}),$ (D.1)

where

$\displaystyle a_{1}(\Phi,v_{1})=\int_{\Omega}\nabla\Phi\cdot\nabla v_{1}dx,$
(D.2) $\displaystyle
a_{2}(\Phi,v_{1})=\frac{1}{\epsilon}\left(\int_{\Gamma_{1}}\Phi
v_{1}ds-\frac{R\sigma}{R\int_{\Gamma_{1}}\sigma
ds+\epsilon}\int_{\Gamma_{1}}\Phi ds\int_{\Gamma_{1}}v_{1}ds\right).$

To see its boundedness property, we have

$\displaystyle|a(\Phi,v_{1})|$
$\displaystyle\leq\|\nabla\Phi\|_{L^{2}(\Omega)}\|\nabla
v_{1}\|_{L^{2}(\Omega)}+\frac{1}{\epsilon}\left(\|\Phi\|_{L^{2}(\Gamma_{1})}\|v_{1}\|_{L^{2}(\Gamma_{1})}+\frac{R\sigma}{R\int_{\Gamma_{1}}\sigma
ds+\epsilon}\|\Phi\|_{L^{1}(\Gamma_{1})}\|v_{1}\|_{L^{1}(\Gamma_{1})}\right)$
(D.3) $\displaystyle\leq\|\nabla\Phi\|_{H^{1}(\Omega)}\|\nabla
v_{1}\|_{H^{1}(\Omega)}+\frac{1}{\epsilon}\left(\|\Phi\|_{L^{2}(\Gamma_{1})}\|v_{1}\|_{L^{2}(\Gamma_{1})}+\frac{R\sigma|\Gamma_{1}|^{2}}{R\int_{\Gamma_{1}}\sigma
ds+\epsilon}\|\Phi\|_{L^{2}(\Gamma_{1})}\|v_{1}\|_{L^{2}(\Gamma_{1})}\right)$
$\displaystyle\leq C_{1}\|\nabla\Phi\|_{H^{1}(\Omega)}\|\nabla
v_{1}\|_{H^{1}(\Omega)},~{}\text{for some constant $C_{1}$ and all
$\Phi,v_{1}\in V_{\Phi}$},$

where the last inequality follows from the trace theorem.

Knowing $\sigma$ is a constant and applying Cauchy Schwarz inequality
$(\int_{\Gamma_{1}}\Phi ds)^{2}\leq|\Gamma_{1}|\int_{\Gamma_{1}}\Phi^{2}ds$,
where $|\Gamma_{1}|$ is the length of $\Gamma_{1}$, we have

$\displaystyle a_{2}(\Phi,\Phi)$
$\displaystyle=\frac{1}{\epsilon}\left(\int_{\Gamma_{1}}\Phi^{2}ds-\frac{R\sigma}{R\int_{\Gamma_{1}}\sigma
ds+\epsilon}\int_{\Gamma_{1}}\Phi ds\int_{\Gamma_{1}}\Phi ds\right)$ (D.4)
$\displaystyle\geq\frac{1}{\epsilon}\bigg{(}1-\frac{R\sigma|\Gamma_{1}|}{R\sigma|\Gamma_{1}|+\epsilon}\bigg{)}\int_{\Gamma_{1}}\Phi^{2}ds$
(D.5) $\displaystyle>0.$

So $a_{2}(\cdot,\cdot)$ is positive definite. By Poincare inequality, we have
that the bilinear form $a_{1}(\cdot,\cdot)$ is strictly positive definite.
Since $a_{1}(\cdot,\cdot)$ is strictly positive definite and
$a_{2}(\cdot,\cdot)$ is positive definite, the bilinear form $a(\cdot,\cdot)$
must be strictly positive definite.

Furthermore, following a similar argument for the bounded property of
$a(\cdot,\cdot)$, we have

$\displaystyle|L(v_{1})|$
$\displaystyle\leq\|f\|_{L^{2}(\Omega)}\|v_{1}\|_{L^{2}(\Omega)}+\frac{1}{\epsilon}\bigg{|}V-\frac{R}{R\int_{\Gamma_{1}}\sigma
ds+\epsilon}\int_{\Gamma_{1}}\sigma Vds\bigg{|}\|v_{1}\|_{L^{1}(\Gamma_{1})}$
(D.6) $\displaystyle\leq C_{2}\|v_{1}\|_{H^{1}(\Omega)},~{}\text{for some
constant $C_{2}$ and all $v_{1}\in V_{\Phi}$}.$

This shows that $L$ is a bounded linear functional on $V_{\Phi}$.

Therefore, the well-posedness follows from the Lax-Milgram lemma [16]. ∎

##  6 Data Availability

The raw data required to reproduce these findings are available to download
from https://data.mendeley.com/ datasets/v6h4tygkvc/1. The processed data
required to reproduce these findings are available to download from
https://data.mendeley.com/datasets/v6h4tygkvc/1.

## References

  * Lee et al. [2008] Y. W. Lee, B.-J. Kim, J.-W. Lim, S. J. Yun, S. Choi, B.-G. Chae, G. Kim, H.-T. Kim,  Metal-insulator transition-induced electrical oscillation in vanadium dioxide thin film,  Applied Physics Letters 92 (2008) 162903.
  * Kim et al. [2010] H.-T. Kim, B.-J. Kim, S. Choi, B.-G. Chae, Y. W. Lee, T. Driscoll, M. M. Qazilbash, D. N. Basov,  Electrical oscillations induced by the metal-insulator transition in vo2,  Journal of Applied Physics 107 (2010) 023702.
  * Shukla et al. [2014] N. Shukla, A. Parihar, E. Freeman, H. Paik, G. Stone, V. Narayanan, H. Wen, Z. Cai, V. Gopalan, R. Engel-Herbert, et al.,  Synchronized charge oscillations in correlated electron systems,  Scientific reports 4 (2014) 4964.
  * Shi and Chen [2018] Y. Shi, L.-Q. Chen,  Phase-field model of insulator-to-metal transition in ${\mathrm{VO}}_{2}$ under an electric field,  Phys. Rev. Materials 2 (2018) 053803.
  * Shi and Chen [2019] Y. Shi, L.-Q. Chen,  Current-driven insulator-to-metal transition in strongly correlated ${\mathrm{vo}}_{2}$,  Phys. Rev. Applied 11 (2019) 014059.
  * Kumar et al. [2017] S. Kumar, J. P. Strachan, R. S. Williams,  Chaotic dynamics in nanoscale $\mathrm{NbO_{2}}$ mott memristors for analogue computing,  Nature 548 (2017) 318–321.
  * Kumar et al. [2020] S. Kumar, R. S. Williams, Z. Wang,  Third-order nanocircuit elements for neuromorphic engineering,  Nature 585 (2020) 518–523.
  * IBabuška [1973] I. IBabuška,  The finite element method with lagrangian multipliers,  Numerische Mathematik 20 (1973) 179–192.
  * Glowinski et al. [1995] R. Glowinski, T.-W. Pan, J. Periaux,  A lagrange multiplier/fictitious domain method for the dirichlet problem—generalization to some flow problems,  Japan journal of industrial and applied mathematics 12 (1995) 87.
  * Nitsche [1971] J. Nitsche,  Ü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 (1971) 9–15.
  * Babuška [1973] I. Babuška,  The finite element method with penalty,  Mathematics of computation 27 (1973) 221–228.
  * ELSGOLC [1961] L. ELSGOLC,  Chapter i - the method of variation in problems with fixed boundaries,  in: L. ELSGOLC (Ed.), Calculus of Variations, International Series of Monographs on Pure and Applied Mathematics, Pergamon, 1961, pp. 13 – 63. URL: http://www.sciencedirect.com/science/article/pii/B9780080095547500066. doi:https://doi.org/10.1016/B978-0-08-009554-7.50006-6.
  * Alnæs et al. [2015] M. S. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, G. N. Wells,  The fenics project version 1.5,  Archive of Numerical Software 3 (2015).
  * Saad and Schultz [1986] Y. Saad, M. H. Schultz,  Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems,  SIAM Journal on Scientific and Statistical Computing 7 (1986) 856–869.
  * Stüben [1999] K. Stüben, Algebraic Multigrid (AMG): An Introduction with Applications, GMD-Forschungszentrum Informationstechnik, 1999.
  * Bressan [2013] A. Bressan, Lecture Notes on Functional Analysis: With Applications to Linear Partial Differential Equations, Graduate studies in mathematics, American Mathematical Society, 2013. URL: https://books.google.com/books?id=dXHQVSSfBnAC.