Solution becomes NaN when parameters are too high

marcbzzz · March 17, 2021, 1:22pm

Hello, I’m quite new to FEniCS and I’m encoutering a problem with my program. I’m calculating the displacement of an elastic membrane under internal pressure.

I can vary different parameters such as the pressure of the number of iterations. When the values are not too big, the program works fine and I get very nice plots.
However if I increase the pressure or the meshSize or the number of iterations too much. The displacement vector becomes NaN everywhere.
And I don’t know exactly why.

Do you have any idea of how I could solve this - as a beginner ?

Here’s my code and the plots I get when the values are too high and when they are the working range

from future import print_function
from fenics import *
import numpy as np
from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
%matplotlib notebook

Constants

Re = 12
Ri = 9
p = 220
dP = 5000
H = 9
Hinv = Re - H

iterations = 89
meshSize = 2
E_tip = 1e5
ratio = 10
E_side = E_tip*ratio
tol = 0.1

Create Domain & Mesh

domain = Circle(Point(0., 0.), Re, 100) - Circle(Point(0., 0.), Ri, 100)
domain = domain - Rectangle(Point(Re, 0.), Point(-Re, -Re))
+ Rectangle(Point(Ri, -10), Point(Re, 0))
+ Rectangle(Point(-Ri,-10),Point(-Re,0))
#domain = domain + Circle(Point(0., 0.), Ri, 100) + R

Create Mesh

mesh = generate_mesh(domain, meshSize)
meshBeforeContraction = generate_mesh(domain, meshSize)
originalMesh = generate_mesh(domain, meshSize)

Set Subdomains and Mark Subdomains

tipDomain = CompiledSubDomain(‘x[0]*x[0] + x[1]x[1] + tol >= RiRi
&& x[0]*x[0] + x[1]x[1] - tol <= ReRe && x[1]>=0- tol &&x[1]>=Hinv-tol
||(x[1] <= 0 + tol &&x[1]>=Hinv && abs(x[0])<=Re+tol && abs(x[0])>=Ri -tol)’,
tol=tol, Hinv=Hinv, Re=Re, Ri=Ri)
sideDomain = CompiledSubDomain('x[1] <= Hinv + tol && x[0]<=Re && x[0]>=Ri -tol ',
tol=tol, Hinv=Hinv, Re=Re, Ri=Ri)

subdomains = MeshFunction(‘size_t’, mesh, mesh.topology().dim(), 0)
subdomains.set_all(4)
tipDomain.mark(subdomains, 5)

Create Boundaries

class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], -10) and on_boundary

class Inner(SubDomain):
def inside(self, x, on_boundary):
return (near(sqrt(x[0]**2+x[1]**2), Ri, 1e-1) and x[1] > Hinv and on_boundary)
or (near(x[0], Ri) and x[1] > -10 and x[1] < 0 and on_boundary)

class InnerHard(SubDomain):
def inside(self, x, on_boundary):
return (near(sqrt(x[0]**2+x[1]**2), Ri, 0.5) and x[1]-tol <= Hinv and x[1] + tol >= 0 and on_boundary)
or (near(abs(x[0]), Ri) and x[1]+tol >= -10 and x[1]-tol <= 0 and x[1]+tol <= Hinv and on_boundary)

class InnerSoft(SubDomain):
def inside(self, x, on_boundary):
return (near(sqrt(x[0]**2+x[1]**2), Ri, 0.5) and x[1] >= Hinv and x[1] >= 0 and on_boundary)
or (near(abs(x[0]), Ri) and x[1] >= -10 and x[1] <= 0 and x[1] >= Hinv and on_boundary)

Mark facets

facets = MeshFunction(“size_t”, mesh, 1)
facets.set_all(0)
Bottom().mark(facets, 1)
Inner().mark(facets, 3)
InnerHard().mark(facets, 31)
InnerSoft().mark(facets, 32)

ds = Measure(“ds”, subdomain_data=facets)

Young Modulus

class Young(UserExpression):
def init(self, subdomains, tip, side, **kwargs):
super().init(**kwargs)
self.subdomains = subdomains
self.side = side
self.tip = tip
def eval_cell(self, values, x, cell):
    if self.subdomains[cell.index] == 4:
        values[0] = self.side
    else:
        values[0] = self.tip
x = SpatialCoordinate(mesh)

Define constants stress/stain :

E = Young(subdomains=subdomains, side=E_side, tip=E_tip, degree=0)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda = Enu/(1+nu)/(1-2nu)

Define Strain Tenosr

def eps(v):
return sym(as_tensor([[v[0].dx(0), 0, v[0].dx(1)],
[0, v[0]/x[0], 0],
[v[1].dx(0), 0, v[1].dx(1)]]))

Define Stress

‘’‘3D hookes law formula of stress’‘’

def sigma(v):
return lmbdatr(eps(v))Identity(3) + 2.0mueps(v)

n = FacetNormal(mesh)

V = VectorFunctionSpace(mesh, ‘CG’, degree=2)
V_original = VectorFunctionSpace(originalMesh, ‘CG’, degree=2)
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*x[0]dx
l = inner(-pn, u_)*x[0]ds(31) + inner(-pn, u_)*x[0]*ds(32)

u = Function(V, name=“Displacement”)
u_original = Function(V_original, name=“original position”)

‘’‘THIS IS THE INTERESTING BOUNDARY CASE’‘’

bcs2 = [DirichletBC(V, Constant((0., 0.)), facets, 1)]

#fid = File(“viz/solution.pvd”)
for i in range(0, iterations):
solve(a == l, u, bcs2)
# u.rename(“u”, “u”) #see the QA reported below.
#fid << u, i
ALE.move(mesh, u)
meshBeforeContraction = Mesh(mesh)
plt.figure()
#plot(mesh, linewidth=0.2)
plot(originalMesh, linewidth=0.2)
p2 = plot(u, mode=“displacement”)
plt.colorbar(p2)
plt.axhline(y=Hinv, color=‘b’, linestyle=‘–’)

plt.show()
plt.title(‘Model before contraction ;’ + 'H = ’ + str(H))

Thank you so much in advance

marcbzzz · March 18, 2021, 4:29pm

Hi, anyone has an idea ? I heard the words “singularity” and preconditioned . However I have no idea of how to implement those here.

dokken · March 19, 2021, 12:17pm

Please format your code appropriately.
You should use 3x` encapsulation around your code, such that quotation, comments and indentation is preserved.

Example:

from dolfin import *
def test():
    mesh = UnitSquareMesh(10, 10)
    V = FunctionSpace(mesh, "CG", 1)
    return V
test()

dparsons · March 19, 2021, 2:19pm

Since you have a working solution in hand, one general approach is to use that solution iteratively. Apply the previous solution as an initial guess for the next iteration, increasing your control parameter incrementally at each iteration until it reaches the value you want.

marcbzzz · March 19, 2021, 4:13pm


from __future__ import print_function
from fenics import *
import numpy as np
from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'notebook')

# Constants
H = 11
p = Constant(100.)
Hinv = 11 - H
Re = 11.
Ri = 9
iterations = 20

¨
# Create Domain & Mesh
domain = Circle(Point(0., 0.), Re, 100) - Circle(Point(0., 0.), Ri, 100)
domain = domain - Rectangle(Point(Re, 0.), Point(-Re, -Re))                 + Rectangle(Point(Ri, -10), Point(Re, 0))                 + Rectangle(Point(-Ri,-10),Point(-Re,0))
# Create Mesh
mesh = generate_mesh(domain, 5)
originalMesh = generate_mesh(domain, 100)

# Set Subdomains and Mark Subdomains

tol = 0.1
tipDomain = CompiledSubDomain('x[0]*x[0] + x[1]*x[1] + tol >= Ri*Ri                                 && x[0]*x[0] + x[1]*x[1] - tol <= Re*Re                                 && x[1]>=0- tol &&x[1]>=Hinv-tol                                  ||(x[1] <= 0 + tol &&x[1]>=Hinv && abs(x[0])<=Re+tol   && abs(x[0])>=Ri -tol)',
                              tol=tol, Hinv=Hinv, Re=Re, Ri=Ri)
sideDomain = CompiledSubDomain('x[1] <= Hinv + tol && abs(x[0])<=Re+tol  && abs(x[0])>=Ri -tol ',
                               tol=tol, Hinv=Hinv, Re=Re, Ri=Ri)


subdomains = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)
subdomains.set_all(4)
tipDomain.mark(subdomains, 5)


# Create Boundaries
class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], -10) and abs(x[0])>=3 and on_boundary


class Outer(SubDomain):
    def inside(self, x, on_boundary):
        return (near(sqrt(x[0]**2+x[1]**2), Ri, 1e-1) and x[1] > Hinv and on_boundary)            or (near(x[0], abs(Ri)) and x[1] > -10 and x[1] < 0 and on_boundary)
    
class TopMiddle(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],0)


# Mark facets
facets = MeshFunction("size_t", mesh, 1)
facets.set_all(0)
Bottom().mark(facets, 1)
Outer().mark(facets, 3)
TopMiddle().mark(facets,2)
ds = Measure("ds", subdomain_data=facets)


# In[3]:


# Plot the mesh and horizontal line
plt.figure()
plot(mesh)

plt.axhline(y=Hinv, color='b', linestyle='--')

plot(subdomains)


class Young(UserExpression):
    def __init__(self, subdomains, tip, side, **kwargs):
        super().__init__(**kwargs)
        self.subdomains = subdomains
        self.side = side
        self.tip = tip

    def eval_cell(self, values, x, cell):
        if self.subdomains[cell.index] == 4:
            values[0] = self.side
        else:
            values[0] = self.tip


x = SpatialCoordinate(mesh)
# Define constants stress/stain :
E_tip = 1e4
E_side = 1e6
E = Young(subdomains=subdomains, side=E_side, tip=E_tip, degree=0)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
tol = 1e-1


# Define Strain Tenosr
def eps(v):
    return sym(as_tensor([[v[0].dx(0), 0, v[0].dx(1)],
                          [0, v[0]/x[0], 0],
                          [v[1].dx(0), 0, v[1].dx(1)]]))


# Define Stress

def sigma(v):
    return lmbda*tr(eps(v))*Identity(3) + 2.0*mu*eps(v)

¨
n = FacetNormal(mesh)


V = VectorFunctionSpace(mesh, 'CG', degree=2)
V_original = VectorFunctionSpace(originalMesh, 'CG', degree=2)
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*x[0]*dx    
l = inner(-p*n, u_)*x[0]*ds(3)

u = Function(V, name="Displacement")
u_original = Function(V_original, name="original position")

# Boundaries and solving
bcs2 = [DirichletBC(V, Constant((0., 0.)), facets, 1)]


fid = File("../Desktop/Viz/solution.pvd")
for i in range(0, iterations):
    solve(a == l, u, bcs2)

    #u.rename("u", "u") #see the QA reported below.
    #fid << u, i
    ALE.move(mesh, u)


plt.figure()
#plot(mesh, linewidth=0.2)
plot(originalMesh, linewidth=0.2)
p2 = plot(u, mode="displacement")
plt.colorbar(p2)
plt.axhline(y=Hinv, color='b', linestyle='--')


plt.show()

marcbzzz · March 19, 2021, 4:13pm

sorry I just fixed it

marcbzzz · March 19, 2021, 4:28pm

What do you mean exactly ?

dparsons · March 19, 2021, 5:02pm

Depending on the specific algorithm, the algorithm can converge more efficiently by commencing from some initial function (a known approximate solution). This is the case in FEniCS, at least with the nonlinear solver, which internally uses Newton iterations to converge to the solution.

The default initial function used in the FEniCS nonlinear solver is the trivial function, zero everywhere. You can override that by assigning a known function to u (the variable u that is provided to the solve(...) function) before running the solve function.

Note that in your example you’re invoking the high level solve function, which is a wrapper around a LinearSolver, not a NonlinearSolver. You can’t assign an initial guess in the case of the LinearSolver. You’d have to rework your code to use the NonlinearSolver.

marcbzzz · March 22, 2021, 1:56pm

I realized that there was a problem with values of x =0, since my plot seems to have diverge giving crazy displacements when x is near 0.

I think there is a division by 0 somewhere. The thing is , I don’t know how to solve it.

Do you have an idea ?

dparsons · March 22, 2021, 6:44pm

You have v[0]/x[0] . Looks like it’s mathematically divergent. Do you expect v/x to be finite at x=0 ?

marcbzzz · March 23, 2021, 9:06pm

Yes I do. However, even replacing x[0] by a number in v[0]/x[0] doesn’t solve the problem. There still is a division by 0 somewhere that I can’t get rid of

dparsons · March 23, 2021, 9:14pm

Your code does not run. Can you make sure it runs cleanly by copy&paste using the python command? Provide a minimal code example. You’ve got references to ipython which are not relevant to the problem and prevent the code from running.

Topic		Replies	Views
How to fix this error(nan error)	1	518	September 22, 2021
Fenics Nan for solving PDEs systems General	2	656	November 17, 2022
Pressure error is too large variational formulation	2	457	June 23, 2020
Problem with inflation of neo-hookean membrane variational formulation	2	646	May 12, 2020
UCM model - Nonlinearvariationalsolver Errors	0	289	June 8, 2021

Solution becomes NaN when parameters are too high

Constants

Create Domain & Mesh

Create Mesh

Set Subdomains and Mark Subdomains

Create Boundaries

Mark facets

Young Modulus

Define constants stress/stain :

Define Strain Tenosr

Define Stress

Related Topics