Coupled linear problem solves to give NaN? (Describes curvature flows)

variational formulation
1

I am trying to solve a coupled linear problem with mixed mixed Dirichlet and Neumann conditions on [0,\infty) \times [0,\pi]. For those interested the equations are used in describing curvature flows on closed genus 0 surfaces. When running the program my solution is immediately NaN. Here is the problem:

(u_1)_t+k'\cot \theta+k=0

(u_2)_t+k''+k=0

where k=u_1-a^2u_2-(b-1) for constants a,b and the dash indicating a derivative w.r.t. \theta.
The boundary conditions are

initial conditions: u_1(0,\theta)=\frac{a^2}{\sqrt{a^2\sin^2\theta+b^2\cos^2\theta}} , u_2(0,\theta)=\frac{a^2b^2}{\left(\sqrt{a^2\sin^2\theta+b^2\cos^2\theta}\right)^3}
at left and right boundary: u_j(t,0)=u_j(t,\pi)=\frac{a^2}{b}, u_j'(t,0)=u_j'(t,\pi)=0, j=1,2.

By application of the chain rule and integration by parts the variational problem can be posed as

\int v_1((u_1)_t+(k_{1}u_1'+k_2u_2')\cot \theta +k)d\theta=0

\int v_2(u_2)_t-v_2'(k_{1}u_1'+k_2u_2')+v_2k d\theta+\int_{\{0\}\cup\{\pi\}}v_2(k_{1}u_1'+k_2u_2') d\theta=0

where k_j is the partial of k w.r.t. u_j.
So after discretising the above problem in the time domain I find that after the first time step my solution is a NaN, I am worried I am messing up the set up of the problem and cannot see how? I have given my code below. It reproduces the error as described.

from fenics import *
import matplotlib.pyplot as plt
import numpy as np

a=5
b=3/2
Intervals=1000
Time=1
dt=Time/Intervals
t=0

#defines the mesh
tol=0.1
upper=np.pi-tol
lower=tol
mesh=IntervalMesh(1000,lower,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],lower))

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

boundaries = MeshFunction('size_t', mesh,0)
#marks the boundaries 
boundaryLeft.mark(boundaries,1)
boundaryRight.mark(boundaries,2)
ds = ds(subdomain_data=boundaries)

#define init and bcs.
u0=Expression(('pow(a,2)/sqrt(pow(a*sin(x[0]),2)+pow(b*cos(x[0]),2))','pow(a*b,2)/pow(sqrt(pow(a*sin(x[0]),2)+pow(b*cos(x[0]),2)),3)'),degree=1,a=a,b=b)
u0=interpolate(u0,V)

#left bc
uL=Expression(('pow(a,2)/b','pow(a,2)/b'),degree=1,t=t,a=a,b=b)

#right bc
uR=Expression(('pow(a,2)/b','pow(a,2)/b'),degree=1,t=t,a=a,b=b)

#neumann boundary conditions
gL=Constant(0)
gL=interpolate(gL,V.sub(0).collapse())

gR=Constant(0)
gR=interpolate(gR,V.sub(0).collapse())

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

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

bc_list=[bcL,bcR]

#now define the variational problem
u=Function(V)
u=u0
u_1,u_2=split(u)
k=u_1-pow(a,2)*u_2-(b-1)
k10=1#derivatives
k01=-pow(a,2)

cot=Expression('pow(tan(x[0]),-1)',degree=1)

u=TrialFunction(V)#make u a trial function for defining bilin form
u_1,u_2=split(u)

F=v_1*(u_1-u0[0]+dt*(k10*u_1.dx(0)+k01*u_2.dx(0))*cot+k*dt)*dx\

+(v_2*u_2-v_2*u0[1]-dt*(k10*u_1.dx(0)+k01*u_2.dx(0))*v_2.dx(0)+dt*k*v_2)*dx\

+dt*(k10*gL.dx(0)+k01*gL.dx(0))*v_2*ds(0)-dt*(k10*gR.dx(0)+k01*gR.dx(0))*v_2*ds(1)
#could I omit this last line since the neumann b.c. gL=gR=0?

r, l = lhs(F), rhs(F) #this is a linear problem

print('r',r) #these do not print as nan
print('l',l)

u=Function(V)#turns u into functions for solver.
u.assign(u0) #assigns an initial guess for u (init condition)
u_1,u_2=split(u)
for n in range(Intervals):
    t+=dt
    gL.t=t
    gR.t=t
    uL.t=t
    uR.t=t
    #solves
    solve(r==l,u,bc_list)
    u0.assign(u)
    print(u0(lower))

As a side note; Is the way I have implemented the neumann boundary condtitions correct? Naively, one might assume that since the Neumann conditions are zero they could be dropped from the variational problem completely. I suspect this is wrong as we have then no way of telling Fenics we wish to implement these conditions, but I could be wrong.

Thanks for anyone willing to have a look at this.

Best
Morgan

2

First of all, as the function 1/tan(x) is singular at at multiple points in your domain, it does not suprise me that you get a nan answer to the PDE. (This means that when cot is interpolated onto a finite element space, one of it values are nan).
This can at least be partially circumvented by letting your mesh not include the zero node, like:

lower=1e-5
upper=2*np.pi
mesh=IntervalMesh(1000,lower,upper)

and evaluate the solution at u0(lower).

3

Hi Jorgen,

You are right, I expect my solution to shrink fast enough to counteract any singular behaviour within the variational problem, had naively thought that this would mean Fenics wouldn’t of had a problem dealing with it - in hindsight I am not surprised either. I have cut the bad parts out of the domain and adjusted the bc’s however it hasn’t helped, I’ve limited my domain to [0.1,\pi-0.1] and no luck. I’ve updated the question to reflect this.

4

I have noticed it is due to my syntax when defining the variational problem.

Syntax was:

F=(...)dx\
 +(...)dx

I thought this implied a line break so not sure why it was causing an issue, however I took it out and it works.

5

I guess the problem is that you have an empty line between the different parts of the form, thus the line break only adds the empty line.

F=F1\
+F2\
+F3

Should work

1 Like