Dolfinx: Time-dependent problem equilibrating instantly

dolfinx
1

Hello,
I am running a 1D reaction-diffusion model on dolfinx; a much simpler version of the old FEniCS example here. With the parameters I have fed it, the system should reach equilibrium on the time scale of a few days. Instead, it is entirely reacted within a fraction of a second. I ran the same problem in an FEA software (FlexPDE) and got the correct results, so it must be not an issue with the problem but with how I have set it up in dolfinx.

Here is the code:

# imports
from mpi4py import MPI
from dolfinx import fem, io, nls, log, mesh, plot, cpp
import numpy as np
import ufl
from petsc4py.PETSc import ScalarType
import matplotlib.pyplot as plt

L = 5e-3 # depth of tray
D = 6.591637557240176e-06 # diffusivity (m2/s)
k = 6e-2 # rate constant (1/s)
C_0 = 0.00133 * 1000 # atmospheric CO2 (mol/m3)
N_0 = 4.018 * 1000 # initial molar density of Ca(OH)2 (mol/m3)
t_final = 60*60*72
t_step = 1

# adjust units of rate constant
k = k / C_0 / 1000 # rate constant 1/(M*s)

domain = mesh.create_interval(MPI.COMM_WORLD, 100, [0, L])
P1 = ufl.FiniteElement('CG', domain.ufl_cell(), 1)
element = ufl.MixedElement([P1, P1])
V = fem.FunctionSpace(domain, element)

# test functions
v1, v2 = ufl.TestFunctions(V)

# fields to solve for
num_subs = V.num_sub_spaces
spaces = []
maps = []
for i in range(num_subs):
    space_i, map_i = V.sub(i).collapse()
    spaces.append(space_i)
    maps.append(map_i)
    
u = fem.Function(V)
u1 = fem.Function(spaces[0])
u2 = fem.Function(spaces[1])

# concentrations from previous timestep
u_n = fem.Function(V)
u_n1, u_n2 = ufl.split(u_n)

# apply timestep variable to each node in the mesh
dt = fem.Constant(domain, ScalarType(t_step))

# turn floats to Constant objects
D_ = fem.Constant(domain, ScalarType(D))
k_ = fem.Constant(domain, ScalarType(k))

# set initial conditions 
u1.vector.array[:] = C_0 * np.ones(u1.vector.array.shape)
u2.vector.array[:] = N_0 * np.ones(u2.vector.array.shape)
u.vector.array[maps[0]] = u1.vector.array
u.vector.array[maps[1]] = u2.vector.array

u1, u2 = ufl.split(u)

# function determining if a node is on the tray top
def on_top_boundary(x):
    return(np.isclose(x[0], L))

# must use collapsed subspace to determine boundary DOFs
V0, submap = V.sub(0).collapse()

# determine boundary DOFs
boundary_dofs = fem.locate_dofs_geometrical((V.sub(0),V0), on_top_boundary)

# apply dirichlet BC to boundary DOFs
bc = fem.dirichletbc(ScalarType(C_0), boundary_dofs[0], V.sub(0))

u.x.array[bc.dof_indices()[0]] = bc.value.value

F = (1/dt)*(u1 - u_n1)*v1*ufl.dx + \
    (1/dt)*(u2 - u_n2)*v2*ufl.dx + \
    D_*ufl.inner(ufl.grad(u1), ufl.grad(v1))*ufl.dx + \
    k_*u1*u2*v1*ufl.dx + \
    k_*u1*u2*v2*ufl.dx

problem = fem.petsc.NonlinearProblem(F, u, bcs=[bc])
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)

# Time-stepping
time_variables = np.zeros((3, int(t_final/t_step)))
t = 0
for n in range(int(t_final/t_step)):
    # Update current time
    t += t_step

    # Solve variational problem for time step
    solver.solve(u)

    # Update previous solution
    u_n.x.array[:] = u.x.array

Here are what the results look like when I plot them up:

Screenshot 2022-04-04 165325
Screenshot 2022-04-04 165325793×783 64.5 KB

Compared to the correct plots:

Screenshot 2022-04-04 165404
Screenshot 2022-04-04 1654041311×985 109 KB

Perhaps there is an issue with how I am applying the time step or updating the solution?

2

Hi,

I think the issue is here:

When you apply these conditions to u, you are setting an initial guess for the solver, not an initial condition. As your code is written currently, the initial condition for Ca(OH)2 is 0; hence why the reaction is complete immediately. To apply an initial condition, you should apply them to u_n, i.e.:

u_n.vector.array[maps[0]] = u1.vector.array
u_n.vector.array[maps[1]] = u2.vector.array
2 Likes