Working Through Poisson-Nernst-Planck equation with time stepping

General
1

Hello dear all,

I’ve been trying to experiement with solving PNP equations with fenics. I have been working with the template in the post " [How to debug the PETSC ERROR?]"

I converted the code to Dolfin because, ultimately I want to be able to solve my own PNP equation and perform optimization with dolfin-adjoint

From my understanding, this template doesn’t perform the time stepping aspect of the solver, can anyone please provide some insight on how to apporach that aspect?

Thank you so much for your assistance.


from dolfin import *
from dolfin_adjoint import *


import numpy as np 


epsilon0  = 8.8541878128e-12 #vacuum permittivity, F/m 
e0        = 1.60217663e-19 #electron charge 
N_A       = 6.0221408e+23 #Avogadro constant. 
kB        = 1.380649e-23 #m^2 kg s^(-2) K^(-1) 
pot       = -1.200000 #electrode potential
T         = 2.980000000000e+02 #temperature
pzc       = 0.000000 #potential of zero charge
potential = (pot - pzc)* e0 / kB / T 
delta_RP  = 3.200000000000e-10 #the thickness of Helmholtz layer
epsilon_RP= 4.170000 #the dielectric constant of Helmholtz layer
epsilon_s = 78.500000 #the dielectric constant of bulk solution

n0        = N_A
D0        = 1.0e-9
lambda_d  = np.sqrt(epsilon_s * epsilon0 * kB * T / e0**2 / n0)
Ld        = lambda_d * 60 #the position of bulk
factor_a  = delta_RP*epsilon_s/epsilon_RP/lambda_d
factor_1  = Ld/lambda_d
factor_2  = Ld*lambda_d/D0
ns_ions   = 2


Length_limit = Ld / lambda_d


# ## splits the mesh into 2048 elements for MPI
domain = IntervalMesh(2048, 0.0, Length_limit)

# ## creates a space of ns_idons + 1 dimensioned piecewise-linear function over the interval mesh 
V = VectorFunctionSpace(domain, "Lagrange", 1, dim=ns_ions+1)



# ## Sets boundary conditions 

tol = 1e-8
right_boundary = lambda x, on_boundary: on_boundary and near(x[0], Length_limit, tol)

bc_H = DirichletBC(V.sub(0), Constant(1.000000000000e-04), right_boundary)
bc_OH = DirichletBC(V.sub(1), Constant(1.000000000000e-04), right_boundary)
bc_phi = DirichletBC(V.sub(2), Constant(0.000000000000e+00), right_boundary)

bcs = [bc_H, bc_OH, bc_phi]


u_n = Function(V)          
u = Function(V)           
u_n_H, u_n_OH, u_n_phi = split(u_n)
u_H, u_OH, u_phi = split(u)
v_H, v_OH, v_phi = TestFunctions(V)




vec_init = Expression(
    ("c0",         
     "c0",         
     "a*x[0] + b"  
    ),
    degree=1,
    c0=1.000000000000e-04,
    a=7.944060176975e-03,
    b=-4.672608459440e+01
)

u_n.interpolate(vec_init)


flux_OH_at_left = Constant(1.0e-5)
dt = Constant(1.0e-3)


flux_H = 9.311000000000e+00*u_H*(1.0)*grad(u_n_phi)+9.311000000000e+00*grad(u_H)
flux_OH = 5.273000000000e+00*u_OH*(-1.0)*grad(u_n_phi)+5.273000000000e+00*grad(u_OH)

F_H = u_H*v_H*dx + dt*factor_1*dot(flux_H,grad(v_H))*dx -u_n_H*v_H*dx
F_OH = u_OH*v_OH*dx + dt*factor_1*dot(flux_OH,grad(v_OH))*dx -dt*factor_1*flux_OH_at_left*v_OH*ds -u_n_OH*v_OH*dx
F_phi = (u_phi - potential)/factor_a*v_phi*ds - dot(grad(u_phi), grad(v_phi))*dx + ((1.0)*u_H+(-1.0)*u_OH)*v_phi*dx

F = F_H + F_OH + F_phi

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

prm = solver.parameters["newton_solver"]
prm["convergence_criterion"] = "residual"  
prm["relative_tolerance"]    = 1e-6    


prm["linear_solver"] = "lu"        

PETScOptions.set("ksp_type", "preonly")  
PETScOptions.set("pc_type", "lu")
PETScOptions.set("pc_factor_mat_solver_type", "mumps") 

solver.solve()