Poisson Nernst Planck with TRBDF2 and SUPG Error

erickmartinez · July 17, 2019, 6:03pm

Hi everyone,

I am trying to solve Poisson-Nernst-Planck equation (PNP). I am attaching a PDF with the formulation for the weak form of the equations.

For the time stepping I implemented a TRBDF-2 rule and to stabilize the velocity I used a SUPG term. Both of these approaches work on the NP problem.

This is the core of the code:

ds  = Measure('ds', domain=mesh, subdomain_data=boundaries)
dS  = Measure('dS', domain=mesh, subdomain_data=boundaries)
dx  = Measure('dx', domain=mesh, subdomain_data=boundaries)
CG1 = FiniteElement("CG", mesh.ufl_cell(), 1)
#Defining the function spaces for the concentration
V_c = FunctionSpace(mesh, 'CG', 1)
#Defining the function spaces for the voltage
V_p = FunctionSpace(mesh, 'CG', 1)

#Defining the mixed function space
W_elem  = MixedElement([CG1, CG1])
W       = FunctionSpace(mesh, W_elem)

#Defining the "Trial" functions
u               = Function(W)      # For time i+1
c, phi          = split(u)         
u_n             = Function(W)      # For time i+1/2
c_n, phi_n      = split(u_n)
u_G             = Function(W)      # For time i
c_G, phi_G      = split(u_G)


v         = TestFunction(W)
(va,vb)   = split(v)

du              = TrialFunction(W)
(dc, dp)        = split(du) 


vde  = Function(V_c)
#tol = 1E-14
Dc   = Expression("D1", D1=D1ums,degree=0)
vel  = Expression("v1", v1=vd1,degree=0)
vde.assign(vel)

qeps = Expression("v1",v1=qe1,degree=0)  


bcs = DirichletBC(W.sub(0),1.0,boundaries,1)#,method="pointwise")

def getVariationalForm(va,vb,p,c):
	a = Dc*(dot(grad(c),grad(va)) - dot(grad(p),grad(vb)) \
		 + zeDkT*c*dot(grad(phi),grad(va)) - zeDkT*Dx(c,0)*Dx(phi,0)*va + vde*Dx(c,0)*va \
		 + qe1*c*vb)*dx
	return a


def getTRBDF2ta(c,phi,i):
	r = Dc*(dot(grad(c),grad(va)) - dot(grad(phi),grad(va)) \
		 + zeDkT*c*dot(grad(phi),grad(va)) - zeDkT*Dx(c,0)*Dx(phi,0)*va + vde*Dx(c,0)*va \
		 + qe1*c*vb)
	return r



hk = CellDiameter(mesh)
# SUPG stabilization
b      = vde-zeDkT*Dx(phi,0)
nb     = sqrt(dot(b,b))
Pek    = nb*hk/(2.0*D1ums)

tau = (hk/(2.0*nb))*((exp(2.0*Pek)+1.0)/(exp(2.0*Pek)-1.0)-1.0/Pek)

GAMMA   = 2.0 - np.sqrt(2)  # 0.59
TRF     = Constant(0.5*GAMMA)
BDF2_T1 = Constant(1.0/(GAMMA*(2.0-GAMMA)))
BDF2_T2 = Constant((1.0-GAMMA)*(1.0-GAMMA)/(GAMMA*(2.0-GAMMA)))
BDF2_T3 = Constant((1.0-GAMMA)/(2.0-GAMMA))

# get the skew symmetric part of the L operator
# LSSNP = dot(vel2,Dx(v2,0))
Lss = (zeDkT*Dx(phi,0) - vde)*Dx(c,0) - (zeDkT/2)*div(grad(phi))
# SUPG Stabilization term
ta = getTRBDF2ta(c_G,phi_G,1); tb = getTRBDF2ta(c_n,phi_n,1)
tc = getTRBDF2ta(c,phi,1)
ra = ((1/dt)*(c_G - c_n) + TRF*( ta + tb) )\
		* tau*Lss*dx
rb  = (c/dt - BDF2_T1*c_G/dt + BDF2_T2*c_n/dt \
	  + BDF2_T3*tc )*tau*Lss*dx

a00 = getVariationalForm(va,vb,phi_n,c_n)
a0L = getVariationalForm(va,vb,phi_G,c_G)
a01 = getVariationalForm(va,vb,phi,c)


A1L = (c_G*va/dt)*dx - (c_n*va/dt)*dx + TRF*(a0L + a00)
A1  = (c/dt - BDF2_T1*c_G/dt + BDF2_T2*c_n/dt)*va*dx \
	+ BDF2_T3*(a01)
F1_G    = A1L + ra
F1      = A1  + rb

J1_G    = derivative(F1_G, u_G, dc)
J1      = derivative(F1, u, dc)
	

ffc_options = {"optimize": True, "quadrature_degree": 8}


newton_solver_parameters = {"nonlinear_solver": "newton",
							"newton_solver" : {
									"linear_solver" : "lu",
									"convergence_criterion": "incremental",
									"absolute_tolerance" : 1E-6,
									"relative_tolerance" : 1E-6,
									"maximum_iterations" : 50,
									"relaxation_parameter" : 1.0 }}

problem1 = NonlinearVariationalProblem(F1, c, bcs, J1, form_compiler_parameters=ffc_options)
problem1G = NonlinearVariationalProblem(F1_G, c_G, bcs1, J1_G, form_compiler_parameters=ffc_options)
solver1  = NonlinearVariationalSolver(problem1)
solver1.parameters.update(newton_solver_parameters)
solver1G = NonlinearVariationalSolver(problem1G)
solver1G.parameters.update(newton_solver_parameters)

t = 0
while t < tmax
	solver1G.solve() 
	solver1.solve()
	(c,phi) = u.split(True)
	c_n.assign(c)
	t += dt

And I am getting the following message

File “./single_layer_PNP.py”, line 512, in update_line
solve(F1_G == 0, u, bcs)
File “/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py”, line 220, in solve
_solve_varproblem(*args, **kwargs)
File “/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py”, line 261, in _solve_varproblem
form_compiler_parameters=form_compiler_parameters)
File “/usr/local/lib/python3.6/dist-packages/dolfin/fem/problem.py”, line 90, in init
F = Form(F, form_compiler_parameters=form_compiler_parameters)
File “/usr/local/lib/python3.6/dist-packages/dolfin/fem/form.py”, line 44, in init
mpi_comm=mesh.mpi_comm())
File “/usr/local/lib/python3.6/dist-packages/dolfin/jit/jit.py”, line 47, in mpi_jit
return local_jit(*args, **kwargs)
File “/usr/local/lib/python3.6/dist-packages/dolfin/jit/jit.py”, line 97, in ffc_jit
return ffc.jit(ufl_form, parameters=p)
File “/usr/local/lib/python3.6/dist-packages/ffc/jitcompiler.py”, line 217, in jit
module = jit_build(ufl_object, module_name, parameters)
File “/usr/local/lib/python3.6/dist-packages/ffc/jitcompiler.py”, line 133, in jit_build
generate=jit_generate)
File “/usr/local/lib/python3.6/dist-packages/dijitso/jit.py”, line 165, in jit
header, source, dependencies = generate(jitable, name, signature, params[“generator”])
File “/usr/local/lib/python3.6/dist-packages/ffc/jitcompiler.py”, line 66, in jit_generate
prefix=module_name, parameters=parameters, jit=True)
File “/usr/local/lib/python3.6/dist-packages/ffc/compiler.py”, line 143, in compile_form
prefix, parameters, jit)
File “/usr/local/lib/python3.6/dist-packages/ffc/compiler.py”, line 185, in compile_ufl_objects
analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
File “/usr/local/lib/python3.6/dist-packages/ffc/analysis.py”, line 90, in analyze_ufl_objects
for form in forms)
File “/usr/local/lib/python3.6/dist-packages/ffc/analysis.py”, line 90, in
for form in forms)
File “/usr/local/lib/python3.6/dist-packages/ffc/analysis.py”, line 174, in _analyze_form
do_apply_restrictions=True)
File “/usr/local/lib/python3.6/dist-packages/ufl/algorithms/compute_form_data.py”, line 418, in compute_form_data
check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode) # Currently testing how fast this is
File “/usr/local/lib/python3.6/dist-packages/ufl/algorithms/check_arities.py”, line 177, in check_form_arity
check_integrand_arity(itg.integrand(), arguments, complex_mode)
File “/usr/local/lib/python3.6/dist-packages/ufl/algorithms/check_arities.py”, line 159, in check_integrand_arity
arg_tuples = map_expr_dag(rules, expr, compress=False)
File “/usr/local/lib/python3.6/dist-packages/ufl/corealg/map_dag.py”, line 37, in map_expr_dag
result, = map_expr_dags(function, [expression], compress=compress)
File “/usr/local/lib/python3.6/dist-packages/ufl/corealg/map_dag.py”, line 86, in map_expr_dags
r = handlers[v.ufl_typecode](v, *[vcache[u] for u in v.ufl_operands])
File “/usr/local/lib/python3.6/dist-packages/ufl/algorithms/check_arities.py”, line 48, in sum
raise ArityMismatch(“Adding expressions with non-matching form arguments {0} vs {1}.”.format(_afmt(a), _afmt(b)))
ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments () vs (‘v_0’,).
Calling FFC just-in-time (JIT) compiler, this may take some time.
Traceback (most recent call last):
File “./single_layer_PNP.py”, line 569, in
main(args.config)
File “./single_layer_PNP.py”, line 553, in main
ani.save(‘EVA’ + ‘.mp4’, writer=writer)
File “/usr/local/lib/python3.6/dist-packages/matplotlib/animation.py”, line 1169, in save
anim._init_draw()
File “/usr/local/lib/python3.6/dist-packages/matplotlib/animation.py”, line 1740, in _init_draw
self._draw_frame(next(self.new_frame_seq()))
File “/usr/local/lib/python3.6/dist-packages/matplotlib/animation.py”, line 1762, in _draw_frame
self._drawn_artists = self._func(framedata, *self._args)
File “./single_layer_PNP.py”, line 512, in update_line
solve(F1_G == 0, u, bcs)
File “/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py”, line 220, in solve
_solve_varproblem(*args, **kwargs)
File “/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py”, line 261, in _solve_varproblem
form_compiler_parameters=form_compiler_parameters)
File “/usr/local/lib/python3.6/dist-packages/dolfin/fem/problem.py”, line 90, in init
F = Form(F, form_compiler_parameters=form_compiler_parameters)
File “/usr/local/lib/python3.6/dist-packages/dolfin/fem/form.py”, line 44, in init
mpi_comm=mesh.mpi_comm())
File “/usr/local/lib/python3.6/dist-packages/dolfin/jit/jit.py”, line 47, in mpi_jit
return local_jit(*args, **kwargs)
File “/usr/local/lib/python3.6/dist-packages/dolfin/jit/jit.py”, line 97, in ffc_jit
return ffc.jit(ufl_form, parameters=p)
File “/usr/local/lib/python3.6/dist-packages/ffc/jitcompiler.py”, line 217, in jit
module = jit_build(ufl_object, module_name, parameters)
File “/usr/local/lib/python3.6/dist-packages/ffc/jitcompiler.py”, line 133, in jit_build
generate=jit_generate)
File “/usr/local/lib/python3.6/dist-packages/dijitso/jit.py”, line 165, in jit
header, source, dependencies = generate(jitable, name, signature, params[“generator”])
File “/usr/local/lib/python3.6/dist-packages/ffc/jitcompiler.py”, line 66, in jit_generate
prefix=module_name, parameters=parameters, jit=True)
File “/usr/local/lib/python3.6/dist-packages/ffc/compiler.py”, line 143, in compile_form
prefix, parameters, jit)
File “/usr/local/lib/python3.6/dist-packages/ffc/compiler.py”, line 185, in compile_ufl_objects
analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
File “/usr/local/lib/python3.6/dist-packages/ffc/analysis.py”, line 90, in analyze_ufl_objects
for form in forms)
File “/usr/local/lib/python3.6/dist-packages/ffc/analysis.py”, line 90, in
for form in forms)
File “/usr/local/lib/python3.6/dist-packages/ffc/analysis.py”, line 174, in _analyze_form
do_apply_restrictions=True)
File “/usr/local/lib/python3.6/dist-packages/ufl/algorithms/compute_form_data.py”, line 418, in compute_form_data
check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode) # Currently testing how fast this is
File “/usr/local/lib/python3.6/dist-packages/ufl/algorithms/check_arities.py”, line 177, in check_form_arity
check_integrand_arity(itg.integrand(), arguments, complex_mode)
File “/usr/local/lib/python3.6/dist-packages/ufl/algorithms/check_arities.py”, line 159, in check_integrand_arity
arg_tuples = map_expr_dag(rules, expr, compress=False)
File “/usr/local/lib/python3.6/dist-packages/ufl/corealg/map_dag.py”, line 37, in map_expr_dag
result, = map_expr_dags(function, [expression], compress=compress)
File “/usr/local/lib/python3.6/dist-packages/ufl/corealg/map_dag.py”, line 86, in map_expr_dags
r = handlers[v.ufl_typecode](v, *[vcache[u] for u in v.ufl_operands])
File “/usr/local/lib/python3.6/dist-packages/ufl/algorithms/check_arities.py”, line 48, in sum
raise ArityMismatch(“Adding expressions with non-matching form arguments {0} vs {1}.”.format(_afmt(a), _afmt(b)))
ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments () vs (‘v_0’,).

If I remove the SUPG term it still gives me the same error. My guess is that there is something wrong with the variational form.

kamensky · July 17, 2019, 6:15pm

The problem is an arity mismatch; it looks like there must be a term in the residual F1_G that does not involve the test function. (The whole residual needs to be linear in the test function for the problem to make sense.) I didn’t try to trace the code in detail, but, taking a quick look, it seems that the term ra does not involve the test function.

erickmartinez · July 17, 2019, 6:26pm

Thank you @kamensky

If I remove the term ‘ra’ and ‘rb’ terms and I got a new trace stack:

Traceback (most recent call last):
File “./single_layer_PNP.py”, line 569, in
main(args.config)
File “./single_layer_PNP.py”, line 366, in main
J1_G = derivative(F1_G, u_G, dc)
File “/usr/local/lib/python3.6/dist-packages/dolfin/fem/formmanipulations.py”, line 82, in derivative
return ufl.derivative(form, u, du, coefficient_derivatives)
File “/usr/local/lib/python3.6/dist-packages/ufl/formoperators.py”, line 281, in derivative
argument)
File “/usr/local/lib/python3.6/dist-packages/ufl/formoperators.py”, line 225, in _handle_derivative_arguments
error(“Coefficient and argument shapes do not match!”)
File “/usr/local/lib/python3.6/dist-packages/ufl/log.py”, line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Coefficient and argument shapes do not match!

So, there seems to be a problem with the derivative.

kamensky · July 17, 2019, 6:33pm

The problem now is that the direction that you gave for the derivative is dc, which is just one component of the trial function, whereas you are differentiating with respect to the full function, u_G (so the “shapes do not match”). Changing dc to du should get rid of this error.

erickmartinez · July 17, 2019, 6:35pm

Thank you @kamensky, it did, but the solver did not converge. I though this had to do with the stability but, when I added the SUPG term its went back to the parity error. I am checking that part.

erickmartinez · July 17, 2019, 10:51pm

As suggested by @kamensky I fixed the direction to match the full function and I managed to fix the terms of the SUPG term.

The problem with the SUPG term is that I had

L_{ss}(c,\phi) , rather than L_{ss}v_a = \left(\frac{zeD}{kT}\nabla\phi - \mathbf{b}\right)\cdot\nabla v_a + \left(\frac{zeD}{2kT}\right)\nabla\cdot\nabla\phi v_a
\left( \left( Lu,v\right) ,\left( L_{ss}(c,\phi )\right) \right) , instead of \left(Lu,L_{ss}v\right)

I fixed those terms and got rid of the warings.

However, I am now stuck with a Newton solver error

Solving nonlinear variational problem.
Newton iteration 0: r (abs) = inf (tol = 1.000e-06) r (rel) = -nan (tol = 1.000e-06)

Topic		Replies	Views
NaN solution problem + general help with my implementation variational formulation	2	770	May 6, 2021
SUPG weak form for a system of reaction-advection-diffusion PDEs	8	1237	January 9, 2022
Solving a Poisson type equation variational formulation	1	471	January 15, 2020
Fokker-Planck solver arity error variational formulation	4	658	February 11, 2020
Non-Linear Coeffecients: dependent on trial functions; variational formulation	2	2648	June 28, 2019

Poisson Nernst Planck with TRBDF2 and SUPG Error

Related Topics