1D time dependent coupled equation: Newton solver fails to converge - what now?

1

Hi all. I am trying to solve the 1D time dependent coupled system u=(u_1,u_2), given by

(u_1)_t=-u_1+u_2(u_2)_xe^{t} and (u_2)_t=\frac{(u_1)_xe^{-t}}{4u_2}-u_2 on (t,x) \in [0,\infty)\times [0,1] which has analytical solution u_1=\sin^2(2x-1)e^{-t}, u_2=\cos(2x-1)e^{-t} (with appropriate Dirichlet boundary conditions).

I have a code which runs, however the newton solver fails to converge even for very small increment in the time step. I am assuming this is happening because of the singular nature of the second equation in the above system. If this failure of convergence is not a result of me incorrectly implementing the program, how am I able to deal with such equations with Fenics? I tried to use Picard iteration but i believe that is already implemented from within the solve() method?

I have included my code below

from fenics import *
from dolfin import*
import matplotlib.pyplot as plt
import numpy as np
#analytical sol: u1=sin^2(2x-t)exp(-t), u2=cos(2x-t)exp(-t)

Intervals=10000
Time=1
dt=Time/Intervals
t=0

#define init and bcs.

#innit
u0=Expression(('pow(sin(2*x[0]),2)','cos(2*x[0])'),degree=1)

#left bc
uL=Expression(('pow(sin(t),2)*exp(-t)','cos(t)*exp(-t)'),degree=1,t=t)

#right bc
uR=Expression(('pow(sin(2-t),2)*exp(-t)','cos(2-t)*exp(-t)'),degree=1,t=t)

#defines the mesh - the unit interval
upper=1
mesh=IntervalMesh(1000,0,upper)
#defines product function space
P1 = FiniteElement('P',interval, 1)
element = MixedElement([P1,P1])
V = FunctionSpace(mesh,element)


#defines boundary classes
class BoundaryRight(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (near(x[0],upper))

class BoundaryLeft(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (near(x[0],0))

#creates instances to use in the program
boundaryLeft = BoundaryLeft()
boundaryRight = BoundaryRight()

#need to create test function as a vector v1,v2
v_1,v_2 = TestFunctions(V)


#this creates the function and then defines its components.
u = Function(V)
u_1,u_2=split(u)


#sets bc at the boundaries
bcL=DirichletBC(V,uL,boundaryLeft)
bcR=DirichletBC(V,uR,boundaryRight)

bc_list=[bcL,bcR]

#now define the variational problem

#interpolates initial cond on the mesh
un=interpolate(u0,V)


F=v_1*(u_1-un[0]+dt*u_1-dt*u_2.dx(0)*u_2*exp(t))*dx\
+v_2*(4*u_2*u_2-4*u_2*un[1]-dt*u_1.dx(0)*exp(-t)+4*u_2*u_2*dt)*dx

for n in range(Intervals):
    t+=dt
    uL.t=t
    uR.t=t
    #solves
    solve(F==0,u,[bcL,bcR])
    un.assign(u)

I am not quite sure what is going on, any advice would be greatly appreciated.
Best
Morgan

2

Maybe consider solving for the variables u_1 and w := u_2^2. Multiplying the second equation through by u_2 (as you did in your code) and noting that

\frac{1}{2}\partial_y w = \frac{1}{2}\partial_y(u_2^2) = u_2(\partial_y u_2)

for y \in \{x,t\}, you’ll see that this actually gets rid of the nonlinearity altogether, since the first equation is

\partial_t u_1 = -u_1 + \frac{1}{2}(\partial_xw)e^t

and the second equation is

\partial_t w = \frac{1}{2}(\partial_xu_1)e^{-t} - 2w\text{ ,}

which is a linear system in the unknowns u_1 and w.

Of course, this first-order system probably also needs some sort of stabilization to get good quality discrete solutions, but that is a separate concern from the nonlinear convergence.

2 Likes
3

Hi Kamensky,

Thanks for your reply! I am quite new to numerics, thanks for showing me this trick - it’s quite clever, and should do the job nicely! If you don’t mind, may I ask you to elaborate on first order systems needing some sort of stabilisation within Fenics?

Thanks so much for your answer.

4

So, your problem fits the general template of multidimensional advection–reaction (where “multidimensional” refers to multiple components of the solution, not multiple space dimensions):

\mathbf{U}_{,t} + \mathbf{A}\cdot\mathbf{U}_{,x} + \mathbf{S}\cdot\mathbf{U} = \mathbf{0}\text{ ,}

where

\mathbf{U} = \left(\begin{array}{c}u_1\\ w\end{array}\right)~~~,~~~\mathbf{A} = -\frac{1}{2}\left(\begin{array}{cc}0&e^t\\ e^{-t} & 0\end{array}\right)~~~\text{, and}~~~\mathbf{S} = \left(\begin{array}{cc}1&0\\ 0 & 2\end{array}\right)\text{ .}

This is a limiting case of advection–diffusion–reaction, where diffusion goes to zero. If diffusion is large relative to advection and reaction, the bilinear form from Galerkin’s method is coercive and bounded in H^1 with similarly-sized constants, so you get a nice H^1 error estimate. However, as diffusion goes to zero, you lose this property, which often manifests as spurious oscillations in discrete solutions.

There are many different finite element methods to handle advection robustly in the literature, and most of the common ones can be implemented in FEniCS. I am partial to the streamline-upwind Petrov–Galerkin (SUPG) method and related formulations, but that may just be an artifact of my academic pedigree. There is an undocumented demo in DOLFIN implementing SUPG for scalar advection–diffusion. I wrote up some notes on the theory of SUPG (following the more recent variational multiscale (VMS) perspective on stabilization) for my graduate class at UCSD. The optimal choice of the SUPG stabilization parameter in multidimensional systems is surprisingly-technical, as derived here, but simpler heuristic choices still work for many systems. Another popular class of methods, especially for problems without diffusion, is the discontinuous Galerkin (DG) approach.

3 Likes
5

Thanks a lot! I very much appreciate thoroughness of this reply. Also thank you for the link to the course notes.