Simulation works in 1D but not in 2D

1

Hi everyone,

I am trying to solve two reaction diffusion equations for two species using mixed elements. I started simulating the problem in 1D and it works fine. When transferring it to 2D, however, the solver solves the first time step and then each other time step gives the same results as the initial step which is not happening in 1D. I am using the 2019 fenics Ubuntu version.

The code below is for the 1D problem. There are two lines that need to be switched to make it 2D (they are provided as comments). The first line is the mesh, the second line is the finite element P1.

If anyone could give me a hint as to what is going on Iā€™d be very thankful.

from dolfin import *
import matplotlib.pyplot as plt

# diffusion and reaction constants 
D1 = 1e-12
D2 = 1e-12
k1 = 2e-1

# time stepping parameters
dt = 1e0
k = Constant(dt)
t_total = 60  # total simulation time
t = 0  # start

# mesh
elements = 100
start = 0
x_end = 1e-4
mesh = IntervalMesh(elements, start, x_end)  # use for 1D
#mesh = RectangleMesh(Point(start,start), Point(x_end,x_end/8), elements, 12)  # use for 2D

# boundary markers
class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], start) and on_boundary


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


sub_domain = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
sub_domain.set_all(0)

left = LeftBoundary()
right = RightBoundary()

left.mark(sub_domain, 0)  # mark left boundary with a 0
right.mark(sub_domain, 1)  # mark right boundary with a 1

dx = Measure("dx", mesh)

# function space
P1 = FiniteElement("P", "interval", 1)  # use for 1D
#P1 = FiniteElement("P", "triangle", 1)  # use for 2D
element = MixedElement([P1, P1])
W = FunctionSpace(mesh, element)
u1 = Function(W)
u1n = Function(W)
v_1, v_2 = TestFunctions(W)
u_1, u_2 = split(u1)

# setting up initial conditions for each species 
i_v = Expression(("x[0]<x_end/2 ? 1:0", "x[0]>x_end/2 ? 1:0"),x_end=x_end,degree=2)
u1n = interpolate(i_v, W)
u_1n, u_2n = split(u1n)

# boundary conditions
bc = DirichletBC(W.sub(1), 1, right)  # c2 = 1, right boundary

# source terms
R = k1 * u_1 * u_2

# weak form
F =  ((u_1 - u_1n) / k) * v_1 * dx + D1 * dot(grad(u_1), grad(v_1)) * dx + R * v_1 * dx  \
      +((u_2 - u_2n) / k) * v_2 * dx + D2 * dot(grad(u_2), grad(v_2)) * dx - R * v_2 * dx 

# solving loop
#plt.figure()
#plt.ion()
concentration = File("output/concentration.pvd")
while t < t_total:
    t += dt
    solve(F==0, u1, bc)
    u1n.assign(u1)
    u_1_, u_2_ = u1.split()
    concentration << (u_1_, t)
    # some line plots
    #plot(u_1, label="liquid concentration")
    #plot(u_2, label="solid concentration")
    #plt.legend()
    #plt.xlabel("x-coordinate (m)")
    #plt.ylabel("concentration")
    #plt.grid(which="both",linestyle='-.')
    #plt.draw()
    #plt.pause(0.001)
2

I think this may be a similar issue to my answer here. Try setting the absolute_tolerance for the Newton solver to 0:

    solve(F==0, u1, bc, solver_parameters={"newton_solver":{"absolute_tolerance":0.0}})
3

Thank you so much! It works fine now.

1 Like