Upadate solutions in time-dependent non linear problem

General
1

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
Faraday = sc.e * sc.Avogadro  # C/mol
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
             * sc.Avogadro  # molecue/mol
             * 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
flux_1 = -(D[0] * ufl.grad(c_1)
           + D[0]/kT * c_1 * z[0] * e * ufl.grad(pot))
flux_2 = -(D[1] * ufl.grad(c_2)
           + 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
Ap = -ufl.inner(ufl.grad(pot), ufl.grad(v_all[0])) * ufl.dx
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?

Thank you in advance,
Matteo

2

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.Functions 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.

3

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))
4

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 :slight_smile:

5

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))

A = -ufl.dot(ufl.grad(c), ufl.grad(v)) * ufl.dx
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.