Hi, everyone. I’m running on Ubuntu2 0.04 LTS (GNU/Linux 4.4.0-19041-Microsoft x86_64) with
Python 3.8.2 (default, Jul 16 2020, 14:00:26)
[GCC 9.3.0] on linux on Windows. I had an isue before but I’ve updated it.
I’m trying to solve the following coupled PDE’s:
\partial u / \partial t = \nabla (\nabla u) + B*(u*\nabla (\nabla h) + \nabla c \nabla h)
-\nabla u * \nabla h - u* \nabla (\nabla h) = \phi1 *u*(h_app - h)
u(t=0) = 1, h(t=0) = 0
\nabla h(x=1) = 0, \nabla h(x = 0) = Bi*h(x=0), \nabla c(x=1) = 0, c(x=0) = 1
u and h are functions of x and t. Bi and h_app constants, B and phi1 0 for x <.1 and 1 for x>=.1
I’ve had issues with the Robin conditions and defining the subdomain parameters, but I think they’re resolved.
For some reason though, when I get a solution to this problem, the solutions are static when they should be time-dependent.
Here’s my code:
Create 1D mesh and define function space
mesh = IntervalMesh(100,0,1)
V = FunctionSpace(mesh,'P',1)
Q = FunctionSpace(mesh,'P',1)
# Define geometry for electrolyte gap and electrode
tol = 1E-14
B_0 = 0
B_1 = 1
bap = Expression('x[0] < 0.1 + tol ? B_0 : B_1', degree=0, \
tol=tol, B_0=B_0, B_1=B_1)
H_0 = 0
H_1 = 1
phi1 = Expression('x[0] < 0.1 + tol ? H_0 : H_1', degree=0, \
tol=tol, H_0=H_0, H_1=H_1)
# Define test and regular functions
h = Function(Q)
w = TestFunction(Q)
u = Function(V)
v = TestFunction(V)
# Define boundary conditions
lyte = 'near(x[0],0)' # Electrode at electrolyte gap
bcu_lyte = DirichletBC(V, Constant(1), lyte) # Changed from 0
# Newman conditions automatically applied at boundaries not specified \
# Derivatives are 0 at the current collector
bcu = [bcu_lyte]
tol = 1E-14
class BoundaryX0(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0, tol)
class BoundaryX1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 1, tol)
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1, 0)
bx0 = BoundaryX0()
bx0.mark(boundary_markers, 0)
bx1 = BoundaryX1()
bx1.mark(boundary_markers, 1)
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
r = Bi
s = 0
g = 0
bch = {0: {'Robin': (r, s)},
1: {'Neumann': 0}}
integrals_N = []
for i in bch:
if 'Neumann' in bch[i]:
if bch[i]['Neumann'] != 0:
g = bch[i]['Neumann']
integrals_N.append(g*w*ds(i))
integrals_R = []
for i in bch:
if 'Robin' in bch[i]:
r, s = bch[i]['Robin']
integrals_R.append(r*(h-s)*w*ds(i))
# Define initial values
h_0 = Constant(0)
h_n = project(h_0, Q)
u_0 = Constant(1)
u_n = project(u_0, V)
# Define expressions used in variational forms
k = Constant(dt)
# Define variational problem for continuity and current
F1 = (u - u_n)*(v/k)*dx + dot(grad(u),grad(v))*dx \
- bap*dot(grad(u), grad(h_n))*v*dx \
- bap*dot( dot(u, div( grad(h_n) )),v ) *dx
F2 = dot( dot(grad(u_n), grad(h)), w )*dx \
+ dot(dot( u_n, div(grad(h)) ), w)*dx \
+ phi1*(u_n/Cb*exp*h_app)*w*dx - phi1*(u_n/Cb*exp*h)*w*dx \
+sum(integrals_R) + sum(integrals_N)
# Create VTK files for visualization output
vtkfile_u = File('twophase/u.pvd')
vtkfile_h = File('twophase/h.pvd')
# Create progress bar
# progress = Progress('Time-stepping')
progress = Progress('Time-stepping', num_steps)
# set_log_level(PROGRESS)
# Time-stepping
t = 0
for n in range(num_steps):
# Update current time
t += dt
# Solve variational problem for time step
solve(F1 == 0, u, bcu)
solve(F2 == 0, h)
# Save solution to file (VTK)
vtkfile_u << (u, t)
vtkfile_h << (h, t)
# Update previous solution
u_n.assign(u)
h_n.assign(h)
# Update progress bar
#progress.update(t / T)
set_log_level(LogLevel.PROGRESS)
progress += 1
set_log_level(LogLevel.ERROR)
I’m wondering if I should define the variational formula differently. I used the _n version of each function in the equation where it’s multipled by the other function’s test function, but I have gotten mixed up on when that’s appropriate.