I am trying to solve the diffusion equation (1D) in two domains, connected by a segregation flux boundary condition:
J = h\left(C_1 - \frac{C_2}{m}\right)
With initial conditions:
C(x,0) = \begin{cases}
1, \quad x \le 10,\\
0, \quad x > 10
\end{cases}
I am using a Crank-Nicolson scheme for time integration:
C^{j+1}_i - C^{j} = \frac{\Delta t}{2} D_i \nabla^2 C^{j+1}_i + \frac{\Delta t}{2} D_i \nabla^2 C^{j}_i \quad \in \Omega_i
and setting up two different meshes, each one with its own function space. I am estimating the concentrations at the boundary points at each integration step to get the segregation flux.
I am setting up the Neumann BC as follows
D_i \int_{\Omega_i} \left(\frac{\partial^2 C_i}{\partial x^2}\right) v_i d\Omega = D_i \int_{\partial\Omega} \left(\frac{C_i}{\partial n}\right) v_i ds - \int_{\Omega_i} \left(\frac{C_i}{\partial x}\right)\left(\frac{v_i}{\partial x}\right) dx
Here is my code:
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
tmax = 5.0 # final time
dt = 1 # The time step
num_steps = round(tmax/dt) +1 # number of time steps
N = 80 # The number of points
dt = tmax/num_steps
D1 = 0.5 # The diffusion coefficient of layer 1
D2 = 0.1 # The diffusion coefficient of layer 2
m = 0.2 # The segregation coefficient
h = 2.0 # The mass transfer coefficient
L1 = 10.0 # The length of layer 1
L2 = 10.0 # The length of layer 2
L = L1 + L2
dx_ = L / N # The x step size
N1 = round(L1 / dx_) + 1
N2 = round(L2 / dx_) + 1
# Create classes for defining parts of the boundaries and the interior
# of the domain
tol = 1E-14
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0, tol) and on_boundary
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], L, tol) and on_boundary
class Interface(SubDomain):
def inside(self, x, on_boundary):
return near(x[0],L1,tol) and on_boundary
# Initialize domains and BC's
top = Top()
bottom = Bottom()
innerBoundary1 = Interface()
innerBoundary2 = Interface()
# Create mesh and define function space
mesh1 = IntervalMesh(N, 0, L1)
mesh2 = IntervalMesh(N, L1, L)
# Initialize mesh function for boundary domains
boundaries1 = MeshFunction("size_t", mesh1, mesh1.topology().dim()-1)
boundaries1.set_all(0)
top.mark(boundaries1,1)
innerBoundary1.mark(boundaries1,2)
boundaries2 = MeshFunction("size_t", mesh2, mesh2.topology().dim()-1)
boundaries2.set_all(0)
innerBoundary2.mark(boundaries2,1)
bottom.mark(boundaries2,2)
ds1 = Measure('ds', domain=mesh1, subdomain_data=boundaries1)
ds2 = Measure('ds', domain=mesh2, subdomain_data=boundaries2)
V1 = FunctionSpace(mesh1, 'P', 1)
V2 = FunctionSpace(mesh2, 'P', 1)
# Define variational problem
u1 = TrialFunction(V1)
v1 = TestFunction(V1)
u1 = Function(V1)
u1_n = Function(V1)
u2 = TrialFunction(V2)
v2 = TestFunction(V2)
u2 = Function(V2)
u2_n = Function(V2)
# Define the discontinuous diffusion coefficient
tol = 1E-14
Dc = Expression("(x[0] <= d + tol) ? D1 : D2",d=L1,D1=D1,D2=D2,tol=tol,degree=0)
# Get the values of the concentrations at the interface
bv1 = (L1)
c1 = u1(bv1)
c2 = u2(bv1)
gval = (h*(c1-c2/m))
g = Constant("%.3e"% (gval))
g_nval = (h*(c1-c2/m))
g_n = Constant("%.3e"% (g_nval))
u1_0 = Expression("cx1",cx1=1.0,degree=0)
u2_0 = Expression("cx2",cx2=0.0,degree=0)
u1 = interpolate(u1_0,V1)
u2 = interpolate(u2_0,V2)
# Get the variational form of Crank-Nicolson integrating scheme
a1 = u1*v1*dx + 0.5*dt*Dc*dot(grad(u1), grad(v1))*dx - 0.5*dt*g*v1*ds1(2)
L1 = u1_n*v1*dx -0.5*dt*Dc*dot(grad(u1_n), grad(v1))*dx + 0.5*dt*g_n*v1*ds1(2)
a2 = u2*v2*dx + 0.5*dt*Dc*dot(grad(u2), grad(v2))*dx - 0.5*dt*g*v2*ds2(1)
L2 = u2_n*v2*dx - 0.5*dt*Dc*dot(grad(u2_n), grad(v2))*dx + 0.5*dt*g_n*v2*ds2(1)
F1 = a1 - L1
F2 = a2 - L2
# Define a figure to plot the result
fig = plt.figure()
fig.set_size_inches(5.5,5.0,forward=True)
ax1 = plt.subplot2grid((1,1), (0,0))
# Time-stepping
t = 0
for n in range(num_steps):
# Update current time
t += dt
print("Solving for time: %.1f s, C1(0) = %.3e, C1(%.2f) = %.3e, C2(%.2f) = %.3e, g=%.3e" %\
(t,u1(0.0),mesh1.coordinates()[len(mesh1.coordinates())-1],c1,\
mesh2.coordinates()[0],c2,gval))
# Save the segregation flux at previous time
g_nval = (h*(c1-c2/m))
g_n = Constant("%.3e"% (g_nval))
# Compute solution
solve(F1 == 0, u1)
solve(F2 == 0, u2)
# Update previous solution
u1_n.assign(u1)
u2_n.assign(u2)
# Get the current segregation flux
c1 = u1(bv1)
c2 = u2(bv1)
gval = (h*(c1-c2/m))
g = Constant("%.3e"% (gval))
# Define the variational problem
a1 = u1*v1*dx + 0.5*dt*Dc*dot(grad(u1), grad(v1))*dx - 0.5*dt*g*v1*ds1(2)
L1 = u1_n*v1*dx -0.5*dt*Dc*dot(grad(u1_n), grad(v1))*dx + 0.5*dt*g_n*v1*ds1(2)
a2 = u2*v2*dx + 0.5*dt*Dc*dot(grad(u2), grad(v2))*dx - 0.5*dt*g*v2*ds2(1)
L2 = u2_n*v2*dx - 0.5*dt*Dc*dot(grad(u2_n), grad(v2))*dx + 0.5*dt*g_n*v2*ds2(1)
F1 = a1 - L1
F2 = a2 - L2
# Plot solution
plot(u1,marker='o')
plot(u2,marker='o')
ax1.set_xlabel("Depth (a.u.)",fontsize=20)
ax1.set_ylabel("Concentration (a.u.)",fontsize=20)
# Hold plot
plt.tight_layout()
plt.show()
fig.savefig('NP_SEGREGATION_FLUX_RESULTS' + '.png' , dpi=300)
Which produces the following output:
Solving for time: 0.8 s, C1(0) = 1.000e+00, C1(10.00) = 0.000e+00, C2(10.00) = 0.000e+00, g=0.000e+00
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.115e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 2.559e-15 (tol = 1.000e-10) r (rel) = 2.296e-15 (tol = 1.000e-09)
Newton solver finished in 1 iterations and 1 linear solver iterations.
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 0.000e+00 (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton solver finished in 0 iterations and 0 linear solver iterations.
Solving for time: 1.7 s, C1(0) = -1.554e-15, C1(10.00) = -4.441e-16, C2(10.00) = 0.000e+00, g=-8.882e-16
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 4.456e-15 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton solver finished in 0 iterations and 0 linear solver iterations.
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 3.701e-16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton solver finished in 0 iterations and 0 linear solver iterations.
Solving for time: 2.5 s, C1(0) = -1.554e-15, C1(10.00) = -4.441e-16, C2(10.00) = 0.000e+00, g=-8.882e-16
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 4.502e-15 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton solver finished in 0 iterations and 0 linear solver iterations.
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 7.402e-16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton solver finished in 0 iterations and 0 linear solver iterations.
Solving for time: 3.3 s, C1(0) = -1.554e-15, C1(10.00) = -4.441e-16, C2(10.00) = 0.000e+00, g=-8.882e-16
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 4.502e-15 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton solver finished in 0 iterations and 0 linear solver iterations.
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 7.402e-16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton solver finished in 0 iterations and 0 linear solver iterations.
Solving for time: 4.2 s, C1(0) = -1.554e-15, C1(10.00) = -4.441e-16, C2(10.00) = 0.000e+00, g=-8.882e-16
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 4.502e-15 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton solver finished in 0 iterations and 0 linear solver iterations.
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 7.402e-16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton solver finished in 0 iterations and 0 linear solver iterations.
Solving for time: 5.0 s, C1(0) = -1.554e-15, C1(10.00) = -4.441e-16, C2(10.00) = 0.000e+00, g=-8.882e-16
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 4.502e-15 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton solver finished in 0 iterations and 0 linear solver iterations.
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 7.402e-16 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton solver finished in 0 iterations and 0 linear solver iterations.
The solution looks like
But this is not the expected behavior of the concentration which should be more like the ones in
But flipped around the y axis.
My questions are:
- Are the Neumann BCs correct for my code?
- I did not set Dirichlet BCs, since this problem was solved for initial values of C on both sides. If I do set C(x=0,t) = 0 I a negative accumulation at the interface, which is not physical. Is FE fine with initial value problems? Other wise, is there something I am missing for the Dirichlet BCs?
- Can this problem be approached with just one mesh and 2 sub meshes? If so, how can the jump at the interface be handled?
- What would be the best choice of Element functions for this problem?
Thanks!
Update:
I noticed that in
it is stated that, for
\begin{split} - \nabla^{2} u &= f \quad {\rm in} \; \Omega, \\ \nabla u \cdot n&= g \quad {\rm on} \ \partial \Omega.\end{split}
\int_{\Omega} f \, {\rm d} x = - \int_{\partial\Omega} g \, {\rm d} s.
is required.
But for this case
\begin{align}-\nabla^2 C^{j+1} &= -\frac{2}{D\Delta t} C^{j+1} + \frac{2}{D\Delta t} C^{j} + \nabla^2 C^{j} \\ \nabla C^{j+1} \cdot n &= \frac{h}{D}\left(C^{j+1}_1 - \frac{C^{j+1}_2}{m}\right)\end{align}
The suggestion here is to use Lagrange mutlipliers. Does it mean I have to split the function space for both meshes?