2D results not similar to 1D for convection diffusion and poisson equation problem

Hello everyone,
I am new to fenics so I need basic help in understanding what points should be kept in mind while changing 1D code to 2D.
The problem I am facing is that given similar inputs, initial conditions and boundary conditions, I cannot extract similar results from a 2D problem.
Below is code for 1D

from __future__ import print_function
from dolfin import *
import dolfin as d
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
import math
import time

from timeit import default_timer as timer
startime = timer()

D_an = Constant(1.0E-7) # m^2/s
D_ca = Constant(1.0E-7) # m^2/s
mu_an = d.Constant(3.9607E-6) # m^2/s*V
mu_ca = d.Constant(3.9607E-6) # m^2/s*V
z_an = Constant(-1.0)
z_ca = Constant(1.0)
z_fc = Constant(-1.0)
Farad = Constant(9.6487E4) # C/mol
eps0 = Constant(8.854E-12) # As/Vm
epsR = Constant(100.0)
Temp = Constant(293.0)
R = Constant(8.3143)

############# mesh part ########################

mesh = IntervalMesh(40000, 0.0, 15.0E-3) # potential value approaches theoretical if nearly 6 million elements used

refinement_cycles = 3
for _ in range(refinement_cycles):
    refine_cells = MeshFunction("bool", mesh, 1)
    refine_cells.set_all(False)
    for cell in cells(mesh):
        if abs(cell.distance(Point(5.0e-3))) < 1.0e-3:
            refine_cells[cell] = True
        elif abs(cell.distance(Point(10.0e-3))) < 1.0e-3:
            refine_cells[cell] = True
            
    mesh = refine(mesh, refine_cells)
    File('mesh_refine.pvd') << mesh
    
def left_boundary(x, on_boundary):
    return on_boundary and near(x[0], 0.0)
def right_boundary(x, on_boundary):
    return on_boundary and near(x[0], 15.0E-3)

################ FE part ########################

P1 = FiniteElement('P', interval, 1)
element = MixedElement([P1, P1, P1])
ME = FunctionSpace(mesh, element)

subdomain = CompiledSubDomain("x[0] >= 5.0E-3 && x[0] <= 10.0E-3")
subdomains = MeshFunction("size_t", mesh, 1)
subdomains.set_all(0)         
subdomain.mark(subdomains, 1)

fc = Constant(2.0) # FCD
V0_r = FunctionSpace(mesh, 'DG', 0)
fc_function = Function(V0_r)
fc_val = [0.0, fc]

help = np.asarray(subdomains.array(), dtype = np.int32)
fc_function.vector()[:] = np.choose(help, fc_val)
#plot(fc_function)
#plt.show()

Sol_c = Constant(1.0)
Poten = 0.0E-3
l_bc_an = DirichletBC(ME.sub(0), Constant(Sol_c), left_boundary)
r_bc_an = DirichletBC(ME.sub(0), Constant(Sol_c), right_boundary)
l_bc_ca = DirichletBC(ME.sub(1), Constant(Sol_c), left_boundary)
r_bc_ca = DirichletBC(ME.sub(1), Constant(Sol_c), right_boundary)
l_bc_psi = DirichletBC(ME.sub(2), Constant(-Poten), left_boundary)
r_bc_psi = DirichletBC(ME.sub(2), Constant(Poten), right_boundary)
bcs = [l_bc_an, r_bc_an, l_bc_ca, r_bc_ca, l_bc_psi, r_bc_psi]

u = Function(ME)

V = FunctionSpace(mesh, P1)

an_int = interpolate(Expression('x[0] >= 5.0E-3 && x[0] <= 10.0E-3 ? 0.0 : 1.0', degree = 1), V)
ca_int = interpolate(Expression('x[0] >= 5.0E-3 && x[0] <= 10.0E-3 ? 2.0 : 1.0', degree = 1), V)
psi_int = interpolate(Expression('x[0] >= 5.0E-3 && x[0] <= 10.0E-3 ? 0.0 : 0.0', degree = 1), V)
assign(u.sub(0), an_int)
assign(u.sub(1), ca_int)
assign(u.sub(2), psi_int)

y0 = np.linspace(0, 0.015, 600)
an_line0 = np.array([u.sub(0)(point) for point in y0])
ca_line0 = np.array([u.sub(1)(point) for point in y0])
psi_line0 = np.array([u.sub(2)(point) for point in y0])

plt.plot(y0, an_line0, 'r:', lw = 1.5)
plt.grid(color='r', linestyle='-', linewidth=0.1)
plt.xlim(0.00, 0.015)
plt.show()
plt.plot(y0, ca_line0, 'r:', lw = 1.5)
plt.grid(color='r', linestyle='-', linewidth=0.1)
plt.xlim(0.00, 0.015)
plt.show()
plt.plot(y0, psi_line0, 'r:', lw = 1.5)
plt.grid(color='r', linestyle='-', linewidth=0.1)
plt.xlim(0.00, 0.015)
plt.show()

an, ca, psi = split(u)
van, vca, vpsi = TestFunctions(ME)

Fan = D_an*(-inner(grad(an), grad(van))*dx - Farad / R / Temp * z_an*an*inner(grad(psi), grad(van))*dx)
Fca = D_ca*(-inner(grad(ca), grad(vca))*dx - Farad / R / Temp * z_ca*ca*inner(grad(psi), grad(vca))*dx)
Fpsi = inner(grad(psi), grad(vpsi))*dx - (Farad/(eps0*epsR))*(z_an*an + z_ca*ca + z_fc*fc_function)*vpsi*dx

F = Fpsi + Fan + Fca

J = derivative(F, u)
problem = NonlinearVariationalProblem(F, u, bcs = bcs, J = J)
solver = NonlinearVariationalSolver(problem)

#print(info(NonlinearVariationalSolver.default_parameters(), True))

#print(solver.default_parameters().items())
solver.parameters["newton_solver"]["linear_solver"] = "mumps"
#solver.parameters["newton_solver"]["maximum_iterations"] = 10
#solver.parameters["newton_solver"]["relative_tolerance"] = 1E-7
#solver.parameters["newton_solver"]["absolute_tolerance"] = 1E-8
#solver.parameters["newton_solver"]["relaxation_parameter"] = 1.0

solver.solve()

an, ca, psi = u.split()

aftersolveT = timer()

y = np.linspace(0, 0.015, 600)
an_line = np.array([an(point) for point in y])
ca_line = np.array([ca(point) for point in y])
psi_line = np.array([psi(point) for point in y])

plt.plot(y, an_line, 'k-', lw = 1.5)
plt.legend(['$c^-$'], loc='best', markerfirst = False, prop={'size': 15})
plt.title("Anion concentration")
plt.xlabel('x[m]')
plt.ylabel('Concentration [mM]')
plt.grid(color='r', linestyle='-', linewidth=0.1)
plt.xlim(0.0, 0.015)
plt.xticks([0.0, 0.005, 0.01, 0.015])
#plt.ylim(0.00, 2.0)
#plt.yticks([0.0, 0.5, 1.0, 1.5, 2.0])
plt.savefig('an_CS', dpi=None, facecolor='w', edgecolor='w', orientation='portrait', papertype=None, format=None, transparent=False, \
            bbox_inches='tight', pad_inches=0.1, frameon=None)
plt.show()

plt.plot(y, ca_line, 'b--', lw = 1.5)
plt.legend(['$c^+$'], loc='best', markerfirst = False, prop={'size': 15})
plt.title("Cation concentration")
plt.xlabel('x[m]')
plt.ylabel('Concentration [mM]')
plt.grid(color='r', linestyle='-', linewidth=0.1)
plt.xlim(0.0, 0.015)
plt.xticks([0.0, 0.005, 0.01, 0.015])
#plt.ylim(0.00, 6.0)
#plt.yticks([0, 1, 2, 3, 4, 5, 6])
plt.savefig('ca_CS', dpi=None, facecolor='w', edgecolor='w', orientation='portrait', papertype=None, format=None, transparent=False, \
            bbox_inches='tight', pad_inches=0.1, frameon=None)
plt.show()
plt.plot(y, psi_line, 'r:', lw = 1.5)
#out = File("psi.pvd")
#out << psi
plt.title("Electric potential")
plt.xlabel('x[m]')
plt.ylabel('Potential [V]')
plt.grid(color='r', linestyle='-', linewidth=0.1)
plt.legend(['$\psi$'], loc='best', markerfirst = False, prop={'size': 15})
plt.xlim(0.0, 0.015)
plt.xticks([0.0, 0.005, 0.01, 0.015])
plt.savefig('psi_CS', dpi=None, facecolor='w', edgecolor='w', orientation='portrait', papertype=None, format=None, transparent=False, \
            bbox_inches='tight', pad_inches=0.1, frameon=None)
plt.show()

print("Number of DOFs: {}".format(ME.dim()))
print("Max. anion concentratrion: {0:1.3f}" .format(max(an_line)))
print("Min. anion concentratrion: {0:1.3f}" .format(min(an_line)))
print("Max. cation concentration: {0:1.3f}" .format(max(ca_line)))
print("Min. cation concentration: {0:1.3f}" .format(min(ca_line)))
print("Max. potential: {0:1.3f}" .format(max(psi_line)))
print("Min. potential: {0:1.3f}" .format(min(psi_line)*1000))

totime = aftersolveT - startime
#print("Start time is : " + str(round(startime, 2)))
#print("After solve time : " + str(round(aftersolveT, 2)))
print("Total time for Simulation : " + str(round(totime)))

And my 2D code is

from __future__ import print_function
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
import dolfin as d

from timeit import default_timer as timer
startime = timer() 

D_an = d.Constant(1.0E-7) # m^2/s
D_ca = d.Constant(1.0E-7) # m^2/s
mu_an = d.Constant(3.9607E-6) # m^2/s*V
mu_ca = d.Constant(3.9607E-6) # m^2/s*V
z_an = d.Constant(-1.0)
z_ca = d.Constant(1.0)
z_fc = d.Constant(-1.0)
Farad = d.Constant(9.6487E4) # C/mol
eps0 = d.Constant(8.854E-12) # As/Vm
epsR = d.Constant(100.0)
Temp = d.Constant(293.0)
R = d.Constant(8.3143)

################################## mesh part ##################################

mesh = d.RectangleMesh(Point(0.0, 0.0), Point(15.0E-3, 15.0E-3), 499, 501)
#File('mesh_orig.pvd') << mesh

refinement_cycles = 1
for _ in range(refinement_cycles):
    refine_cells = d.MeshFunction("bool", mesh, 2)
    refine_cells.set_all(False)
    for cell in cells(mesh):
        mp = cell.midpoint()
        if mp.x() > 3.0e-3 and mp.x() < 12.0e-3:
            if abs(mp.y() - 10.0e-3) < 2.0E-3:
                refine_cells[cell] = True
            elif abs(mp.y() - 5.0e-3) < 2.0E-3:
                refine_cells[cell] = True
        if mp.y() > 3.0e-3 and mp.y() < 12.0e-3:
            if abs(mp.x() - 10.0e-3) < 2.0E-3:
                refine_cells[cell] = True
            elif abs(mp.x() - 5.0e-3) < 2.0E-3:
                refine_cells[cell] = True
                
    mesh = refine(mesh, refine_cells)
    File('mesh_refine.pvd') << mesh

def left_boundary(x, on_boundary):
    return on_boundary and near(x[0], 0.0)
def right_boundary(x, on_boundary):
    return on_boundary and near(x[0], 15.0E-3)
def bottom_boundary(x, on_boundary):
    return on_boundary and near(x[1], 0.0)
def top_boundary(x, on_boundary):
    return on_boundary and near(x[1], 15.0E-3) 

################################### FE part ###################################   
    
P1 = FiniteElement('P', 'triangle', 2)
element = MixedElement([P1, P1, P1])
ME = FunctionSpace(mesh, element)

subdomain = CompiledSubDomain("x[0] >= 5.0E-3 && x[0] <= 10.0E-3 && x[1] >= 5.0E-3 && x[1] <= 10.0E-3")
subdomains = MeshFunction("size_t", mesh, 2)
subdomains.set_all(0)        
subdomain.mark(subdomains, 1)

fc = d.Constant(2.0) # FCD
V0_r = FunctionSpace(mesh, 'DG', 0)
fc_function = Function(V0_r)
fc_val = [0.0, fc]

help = np.asarray(subdomains.array(), dtype = np.int32)
fc_function.vector()[:] = np.choose(help, fc_val)
zeroth = plot(fc_function, title = "Fixed charge density, $c^f$")
plt.colorbar(zeroth)
plot(fc_function)
plt.xlim(0.0, 0.015)
plt.xticks([0.0, 0.005, 0.01, 0.015])
plt.ylim(0.0, 0.015)
plt.yticks([0.0, 0.005, 0.01, 0.015])
plt.show()

Sol_c = d.Constant(1.0)
Poten = 0E-3
l_bc_an = d.DirichletBC(ME.sub(0), Sol_c, left_boundary)
r_bc_an = d.DirichletBC(ME.sub(0), Sol_c, right_boundary)
l_bc_ca = d.DirichletBC(ME.sub(1), Sol_c, left_boundary)
r_bc_ca = d.DirichletBC(ME.sub(1), Sol_c, right_boundary)
l_bc_psi = DirichletBC(ME.sub(2), Constant(-Poten), left_boundary)
r_bc_psi = DirichletBC(ME.sub(2), Constant(Poten), right_boundary)
bcs = [l_bc_an, r_bc_an, l_bc_ca, r_bc_ca, l_bc_psi, r_bc_psi]

u = Function(ME)

#V = FunctionSpace(mesh, P1)
#an_int = interpolate(Expression('x[0] >= 5.0E-3 and x[0] <= 10.0E-3 && x[1] >= 5.0E-3 and x[1] <= 10.0E-3 ? 0.4142 : 1.0', degree = 1), V)
#ca_int = interpolate(Expression('x[0] >= 5.0E-3 and x[0] <= 10.0E-3 && x[1] >= 5.0E-3 and x[1] <= 10.0E-3 ? 2.4142 : 1.0', degree = 1), V)
#psi_int = interpolate(Expression('x[0] >= 5.0E-3 and x[0] <= 10.0E-3 && x[1] >= 5.0E-3 and x[1] <= 10.0E-3 ? -22.252E-3 : 0.0', degree = 1), V)
#assign(u.sub(0), an_int)
#assign(u.sub(1), ca_int)
#assign(u.sub(2), psi_int)

class InitialConditions(UserExpression):
 def eval(self, values, x):
     if((x[0] >= 5.0E-3) and (x[0] <= 10.0E-3) and (x[1] >= 5.0E-3) and (x[1] <= 10.0E-3)):
         values[0] = 0.0
         values[1] = 2.0
         values[2] = 0.0
#         print("counter value",values)
     else:
         values[0] = 1.0
         values[1] = 1.0
         values[2] = 0.0
#         print("counter value 2 ",values)
 def value_shape(self):
   return (3,)
u.interpolate(InitialConditions(degree=1))

an, ca, psi = split(u)
van, vca, vpsi = TestFunctions(ME)

Fan = D_an*(-inner(grad(an), grad(van))*dx - Farad / R / Temp * z_an*an*inner(grad(psi), grad(van))*dx)
Fca = D_ca*(-inner(grad(ca), grad(vca))*dx - Farad / R / Temp * z_ca*ca*inner(grad(psi), grad(vca))*dx)
Fpsi = inner(grad(psi), grad(vpsi))*dx - (Farad/(eps0*epsR))*(z_an*an + z_ca*ca + z_fc*fc_function)*vpsi*dx

F = Fpsi + Fan + Fca

J = derivative(F, u)
problem = NonlinearVariationalProblem(F, u, bcs = bcs, J = J)
solver = NonlinearVariationalSolver(problem)
#solver = AdaptiveNonlinearVariationalSolver(problem, M)

#print(solver.default_parameters().items())
solver.parameters["newton_solver"]["linear_solver"] = "mumps"
#solver.parameters["newton_solver"]["maximum_iterations"] = 10
#solver.parameters["newton_solver"]["relative_tolerance"] = 1E-7
#solver.parameters["newton_solver"]["absolute_tolerance"] = 1E-5
#solver.parameters["newton_solver"]["relaxation_parameter"] = 1.0
solver.solve()

#########################  Post - Processing  #################################
an, ca, psi = u.split()

aftersolveT = timer() 

first = plot(an, title = "Anion concentration, $c^-$")
plt.colorbar(first)
plt.xlim(0.0, 0.015)
plt.xticks([0.0, 0.005, 0.01, 0.015])
plt.ylim(0.0, 0.015)
plt.yticks([0.0, 0.005, 0.01, 0.015])
plt.show()

second = plot(ca, title = "Cation concentration, $c^+$")
plt.xlim(0.0, 0.015)
plt.xticks([0.0, 0.005, 0.01, 0.015])
plt.ylim(0.0, 0.015)
plt.yticks([0.0, 0.005, 0.01, 0.015])
plt.colorbar(second)
plt.show()

third = plot(psi, title = "Electric potential, $\psi$")
plt.xlim(0.0, 0.015)
plt.xticks([0.0, 0.005, 0.01, 0.015])
plt.ylim(0.0, 0.015)
plt.yticks([0.0, 0.005, 0.01, 0.015])
plt.colorbar(third)
plt.show()

d.File('an_CS.pvd') << an
d.File('ca_CS.pvd') << ca
d.File('psi_CS.pvd') << psi
y_position = 7.5E-3
x_min = 0.0
x_max = 15.0E-3
N_steps = 300
delta_x = (x_max - x_min)/float(N_steps) 

an_vector = []
for i in range(N_steps):
    x_position = float(i * delta_x)
    an_vector.append(an(x_position, y_position))
x_values = [i * delta_x for i in range(N_steps)]
plt.plot(x_values, an_vector, color = 'k', linestyle = '--', linewidth = 1.0)
plt.title("Anion concentration")
plt.xlim(0.0, 0.015)
plt.xticks([0.0, 0.005, 0.01, 0.015])
#plt.ylim(0.0, 2.0)
#plt.yticks([0.0, 0.5, 1.0, 1.5, 2.0])
plt.grid(color = 'k', linestyle = '-', linewidth = 0.1)
plt.legend(['$c^-$'], loc='best', markerfirst = False, prop={'size': 15})
plt.savefig('an_CS_1D')
plt.show()
ca_vector = []
for i in range(N_steps):
    x_position = float(i * delta_x)
    ca_vector.append(ca(x_position, y_position))
x_values = [i * delta_x for i in range(N_steps)]
plt.plot(x_values, ca_vector, color = 'k', linestyle = '-', linewidth = 1.0)
plt.xlim(0.0, 0.015)
plt.xticks([0.0, 0.005, 0.01, 0.015])
plt.title("Cation concentration")
plt.legend(['$c^+$'], loc='best', markerfirst = False, prop={'size': 15})
plt.grid(color = 'k', linestyle = '-', linewidth = 0.1)
plt.savefig('ca_CS_1D')
plt.show()

psi_vector = []
for i in range(N_steps):
    x_position = float(i * delta_x)
    psi_vector.append(psi(x_position, y_position))
x_values = [i * delta_x for i in range(N_steps)]
plt.plot(x_values, psi_vector, color = 'k', linestyle = '-', linewidth = 1.0)
#plt.plot(x_values, psi_vector)
plt.xlim(0.0, 0.015)
plt.xticks([0.0, 0.005, 0.01, 0.015])
plt.grid(color = 'k', linestyle = '-', linewidth = 0.1)
plt.legend(['$\psi$'], loc='best', markerfirst = False, prop={'size': 15})
plt.savefig('psi_CS_1D')
plt.show()

print("Number of DOFs: {}".format(ME.dim()))
print("Max. anion concentratrion: {0:1.3f}" .format(max(an_vector)))
print("Min. anion concentratrion: {0:1.3f}" .format(min(an_vector)))
print("Max. cation concentration: {0:1.3f}" .format(max(ca_vector)))
print("Min. cation concentration: {0:1.3f}" .format(min(ca_vector)))
print("Max. potential: {0:1.3f}" .format(max(psi_vector)))
print("Min. potential: {0:1.3f}" .format(min(psi_vector)))
totime = aftersolveT - startime
print("Total time for Simulation : " + str(round(totime)))
1 Like