# Upadate solutions in time-dependent non linear problem

Dear All,
I’m trying to solve Poisson-Nernst-Planck equation but I’ve problem updating the previous solutions at each time step.

Here my code:

``````from mpi4py import MPI
from petsc4py import PETSc

from dolfinx import default_scalar_type as dst
from dolfinx import mesh, fem, io, log
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver

import ufl
import basix

import scipy.constants as sc
import numpy as np
import pathlib

cache_dir = pathlib.Path('cache_PNP')

# Define quantities
pot_0 = 25 * 1e-3  # V
c_bulk = [0.1 * 1e3, 0.1 * 1e3]  # mol/m^3
z = [1, -1]  # charge/molecule
# from https://www.aqion.de/site/diffusion-coefficients
# NaCl
D = [1.330 * 1e-9, 2.030 * 1e-9]  # m^2/s
kT = sc.Boltzmann * 298.15  # J
eps = sc.epsilon_0 * 78.36  # C/V*m

T = 10  # s
num_step = 100
dt = T/num_step  # s

Is = 0
for c_i, z_i in zip(c_bulk, z):
Is += 1/2 * c_i * z_i**2  # mol*charge / molecule*m^3

k = ufl.sqrt(sc.e**2  # C^2/charge
* 2 * Is  # mol*charge / molecule*m^3
/ (eps * kT)  # V*m/C * 1/J
)  # 1/m

d_l = 1/k  # m

# Create a simple 1D domain
L = 20 * d_l  # m
domain = mesh.create_interval(MPI.COMM_WORLD, 50, [0.0, L])

# Mixed elements
P2 = basix.ufl.element('Lagrange', 'interval', 2)
elements = [P2, P2, P2]  # pot, c_1, c_2
elements_n = [P2, P2]  # c_1_n, c_2_n

element = basix.ufl.mixed_element(elements)
element_n = basix.ufl.mixed_element(elements_n)

# Function space
V = fem.functionspace(domain, element)
W = fem.functionspace(domain, element_n)

# Find boundary factes and dofs on V.sub(0)
tdim = domain.topology.dim
fdim = tdim - 1

def left_boundary(x):
return np.isclose(x[0], 0.0)

left_boundary_facets = mesh.locate_entities_boundary(domain, fdim,
left_boundary)
left_boundary_dofs_0 = fem.locate_dofs_topological(V.sub(0), fdim,
left_boundary_facets)

# Create trial and test functions
u_all = fem.Function(V)  # current solution
v_all = ufl.TestFunction(V)

u_all_list = ufl.split(u_all)
pot = u_all_list[0]
c_1 = u_all_list[1]
c_2 = u_all_list[2]

u_all_n = fem.Function(W)  # for previous step
u_all_n_list = u_all_n.split()
c_1_n = u_all_n_list[0]
c_2_n = u_all_n_list[1]

# Constant declaration
eps = fem.Constant(domain, dst(eps))
kT = fem.Constant(domain, dst(kT))
e = fem.Constant(domain, dst(sc.e))
dt = fem.Constant(domain, dst(dt))
for i, c_i in enumerate(c_bulk):
c_bulk[i] = fem.Constant(domain, dst(c_i))
for i, z_i in enumerate(z):
z[i] = fem.Constant(domain, dst(z_i))
for i, D_i in enumerate(D):
D[i] = fem.Constant(domain, dst(D_i))

# Define flux
+ D[0]/kT * c_1 * z[0] * e * ufl.grad(pot))
+ D[1]/kT * c_2 * z[1] * e * ufl.grad(pot))

# Define source term
rho = Faraday * (z[0] * c_1 + z[1] * c_2)
f = -rho / eps

# Define initial condition
# concentration equal to the bulk conc. everywhere at t=0
c_1_0_expr = fem.Expression(c_bulk[0],
W.sub(0).element.interpolation_points())
c_2_0_expr = fem.Expression(c_bulk[1],
W.sub(1).element.interpolation_points())
# pot_n.interpolate(pot_0_expr)
c_1_n.interpolate(c_1_0_expr)
c_2_n.interpolate(c_2_0_expr)

# ===Weak form===
# Poisson
lp = f * v_all[0] * ufl.dx

# Continuity
Ac_1 = ufl.inner(flux_1, ufl.grad(v_all[1])) * ufl.dx
lc_1 = ((c_1 - c_1_n) / dt) * v_all[1] * ufl.dx
Ac_2 = ufl.inner(flux_2, ufl.grad(v_all[2])) * ufl.dx
lc_2 = ((c_2 - c_2_n) / dt) * v_all[2] * ufl.dx

# BC
bc_0 = fem.dirichletbc(dst(pot_0), left_boundary_dofs_0,
V.sub(0))

# Total probelm
F = (Ap - lp
+ Ac_1 - lc_1
+ Ac_2 - lc_2)

# Solver
problem = NonlinearProblem(F, u_all, [bc_0],
jit_options={'cache_dir': cache_dir})
solver = NewtonSolver(MPI.COMM_WORLD, problem)

solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.max_it = 100
solver.report = True

log.set_log_level(log.LogLevel.INFO)

# Time loop
t = 0

# P1 element for XDMF file
V1 = fem.functionspace(domain, ('Lagrange', 1))
pot_P1 = fem.Function(V1)
c_1_P1 = fem.Function(V1)
c_2_P1 = fem.Function(V1)

pot_out = io.XDMFFile(domain.comm,
'pot.xdmf',
'w')
c_1_out = io.XDMFFile(domain.comm,
'c_1.xdmf',
'w')
c_2_out = io.XDMFFile(domain.comm,
'c_2.xdmf',
'w')
pot_out.write_mesh(V1.mesh)
c_1_out.write_mesh(V1.mesh)
c_2_out.write_mesh(V1.mesh)

for step in range(num_step):
t = t + dt
print(f'Time: {t} (s)')
n, converged = solver.solve(u_all)
if converged:
u_list = u_all.split()
c_1_n.interpolate(u_list[1])
c_2_n.interpolate(u_list[2])

# Save
pot_P1_ex = fem.Expression(u_list[0],
V1.element.interpolation_points())
pot_P1.interpolate(pot_P1_ex)
pot_out.write_function(
pot_P1,
t=t)
c_1_P1_ex = fem.Expression(u_list[1],
V1.element.interpolation_points())
c_1_P1.interpolate(c_1_P1_ex)
c_1_out.write_function(
c_1_P1,
t=t)
c_2_P1_ex = fem.Expression(u_list[2],
V1.element.interpolation_points())
c_2_P1.interpolate(c_2_P1_ex)
c_2_out.write_function(
c_2_P1,
t=t)
else:
print('Not converge')
break

pot_out.close()
c_1_out.close()
c_2_out.close()
``````

When I check the results, I notice that I have obtained constant functions, and these do not vary over time. I believe the problem is that the concentrations at the previous time, c_1_n and c_2_n, are not being updated within the loop. Could someone kindly help me?

Matteo

I have not run the code and so that may not be the actual issue, but I would try defining `W` as a P2 space, rather than a mixed [P2, P2].
In turn, `c_1_n` and `c_2_n` would be `fem.Function`s on this (new) `W`.

Regardless of whether this fixes your issue or not, I’d still suggest you stick with the `W` I propose. Using a mixed space is unnecessary for `u_all_n`: in contrast to `u_all` (which you use as unknown in your problem) you are never solving for `u_all_n`, so it won’t make any difference if `c_1_n` and `c_2_n` are grouped together (in a mixed function space) or not.

This does not resolve the problem of non-evolving solution, but I would interpolate the subspace of the solution into `c_1_n` and `c_2_n` as follows:

``````
c_1_n = fem.Function(fem.functionspace(domain, P2))
c_2_n = fem.Function(fem.functionspace(domain, P2))

for step in range(num_step):
t = t + dt
print(f'Time: {t} (s)')
n, converged = solver.solve(u_all)
if converged:
u_list = u_all.split()
#c_1_n.interpolate(u_list[1])
#c_2_n.interpolate(u_list[2])

c_1_n.interpolate(u_all.sub(1))
c_2_n.interpolate(u_all.sub(2))
``````

Thanks @francesco-ballarin for your answer, I tried to follow your advice but the solution still doesn’t evolving in time. But you’re right, so I’ll keep the change you proposed in my code

Thanks, Jan. The issue might not be related to interpolation. To test a simpler scenario, I wrote a code to solve Fick’s second law:

``````from mpi4py import MPI

from dolfinx import default_scalar_type as dst
from dolfinx import mesh, fem, io, log
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver

import ufl
import basix

import numpy as np
import pathlib

cache_dir = pathlib.Path('cache_Fick')

T = 150
num_step = 200
dt = T/num_step

c_0 = 1.0
c_L = 0.0
c_init = 0.0

L = 20

# Create domain and element
domain = mesh.create_interval(MPI.COMM_WORLD, 50, [0.0, L])
P1 = basix.ufl.element('Lagrange', 'interval', 1)

# Function space
V = fem.functionspace(domain, P1)

# Get boundary dofs
tdim = domain.topology.dim
fdim = tdim - 1

def left_boundary(x):
return np.isclose(x[0], 0.0)

left_boundary_facets = mesh.locate_entities_boundary(domain, fdim,
left_boundary)
left_boundary_dofs = fem.locate_dofs_topological(V, fdim,
left_boundary_facets)

def right_boundary(x):
return np.isclose(x[0], L)

right_boundary_facets = mesh.locate_entities_boundary(domain, fdim,
right_boundary)
right_boundary_dofs = fem.locate_dofs_topological(V, fdim,
right_boundary_facets)

# Create functions
c = fem.Function(V)  # current solution
v = ufl.TestFunction(V)

c_n = fem.Function(V)  # previous solution

# Initial condition
c_init_expr = fem.Expression(fem.Constant(domain,
dst(c_init)),
V.element.interpolation_points(),
jit_options={'cache_dir': cache_dir})
c_n.interpolate(c_init_expr)
c_n.x.scatter_forward()

# Strong form: D * nabla^2 c = dc/dt
# Weak form:
# -D int_omega (nabla c, nabla v) domega = int_omega dc/dt * v domega
# time discretization: dc/dt = (c - c_n) / dt
# assume D = 1
dt = fem.Constant(domain, dst(dt))

lf = ((c - c_n) / dt) * v * ufl.dx

# BC
bc_0 = fem.dirichletbc(dst(c_0), left_boundary_dofs, V)
bc_L = fem.dirichletbc(dst(c_L), right_boundary_dofs, V)

# Problem
F = A - lf

problem = NonlinearProblem(F, c, [bc_0, bc_L],
jit_options={'cache_dir': cache_dir})

# Solver
solver = NewtonSolver(MPI.COMM_WORLD, problem)

solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.max_it = 100
solver.report = True

# Output
c_out = io.XDMFFile(domain.comm,
'fick_out.xdmf',
'w')
c_out.write_mesh(domain)

# Time loop
t = 0
while t < T:
if t == 0:
c_out.write_function(c, t)

t += dst(dt)
print(f'Time {t}')
n, converged = solver.solve(c)
print(f'converged: {converged}, n_iter: {n}')
c_n.interpolate(c)  # update prev. sol
c_n.x.scatter_forward()

# Save
c_out.write_function(c, t)

c_out.close()
``````

In this case, everything works, and the solution evolves over time. Now I belive that the issue with my code for the PNP equation might be related to the weak form or the incorrect use of mixed element spaces.