When I tried to write the density variable obtained by the conditional expression calculation into the pvd file, I got an error message. Besides, it is feasible to draw it by plot command. Can anyone solve this problem?
This is my code:
import random
from numpy.ma import tanh
from matplotlib import pyplot as plt
from ufl import conditional
from fenics import *
import numpy as np
def left(c_):
return (rho_eq_ox - rho_int_ox) / c_ox * c_ + rho_int_ox
def mid(c_):
return (rho_eq_met - rho_eq_ox) / (c_met - c_ox) * (c_ - c_ox) + rho_eq_ox
def right(c_):
return (rho_eq_met - rho_int_met) / (c_met - 1) * (c_ - c_met) + rho_eq_met
class InitialConditions(UserExpression):
def eval(self, values, x):
x_ = x[0]
y_ = x[1]
random.seed(x_)
rx = 2 * random.random() - 1
values[0] = 0.5 * (1 + tanh(2 / wint * (y_ - L - 0.5e-4 * L * rx)))
values[1] = 0.0
def value_shape(self):
return 2,
# Sub domain for Periodic boundary condition
class PeriodicBoundaryX(SubDomain):
def inside(self, x, on_boundary):
return DOLFIN_EPS > x[0] > (- DOLFIN_EPS) and on_boundary
def map(self, x, y):
y[0] = x[0] - L
y[1] = x[1]
# Class for interfacing with the Newton solver
class CahnHilliardEquation(NonlinearProblem):
def __init__(self, a, L):
NonlinearProblem.__init__(self)
self.L = L
self.a = a
def F(self, b, x):
assemble(self.L, tensor=b)
def J(self, A, x):
assemble(self.a, tensor=A)
#
L = 0.176
wint = L / 16
rho_int_ox = 8.24783e3
rho_int_met = 6.96642e3
rho_eq_ox = 7.97994e3
rho_eq_met = 9.01478e3
c_ox = 1.87512e-3
c_met = 6.16085e-1
nx = 192
ny = 384
Vm = 1e-5
R = 8.31446
T = 3000
w_ref = 1e-9
DU = 1.55E-12
DO = 2.99e-9
DZ = DU
DF = DU
M_ref = Vm / (R * T) * (DU + DZ + DO + DF)
M = wint / w_ref * M_ref
gama = 0.3854
A0 = 12 * gama / (c_met - c_ox) ** 4 / wint
dt = 0.02
Episino = 3 * wint * gama / 2 / (c_met - c_ox) ** 2
coef = Episino
mesh = RectangleMesh(Point(0, 0), Point(L, 2 * L), nx, ny)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1 * P1)
du = TrialFunction(ME)
q, v = TestFunctions(ME)
u = Function(ME) # current solution
u0 = Function(ME) # solution from previous converged step
dc, dmu = split(du)
c, mu = split(u)
c0, mu0 = split(u0)
u_init = InitialConditions(degree=1)
u.interpolate(u_init)
u0.interpolate(u_init)
# Compute the chemical potential df/dc
c = variable(c)
f = A0 * (c - c_ox) ** 2 * (c - c_met) ** 2
df_dc = diff(f, c)
L0 = c * q * dx - c0 * q * dx + dt * M * dot(grad(mu), grad(q)) * dx
L1 = mu * v * dx - df_dc * v * dx - coef * dot(grad(c), grad(v)) * dx
L = L0 + L1
a = derivative(L, u, du)
problem = CahnHilliardEquation(a, L)
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6
# Output file
pvd_u = File("u.pvd", "compressed")
pvd_rho = File("rho.pvd", "compressed")
t = 0.0
solver.solve(problem, u.vector())
pvd_u << u.split()[0]
u_0 = u.split(deepcopy=True)[0]
rho = conditional(lt(u_0, c_ox), left(u_0), conditional(lt(u_0, c_met), mid(u_0), right(u_0)))
#pvd_rho << rho #This line is the wrong place
output_rho0 = plot(rho, mode='color')
plt.colorbar(output_rho0)
plt.title("Density")
plt.show()
The error message is:
Traceback (most recent call last):
File "/home/wangyao/PycharmProjects/hello_word/main.py", line 156, in <module>
pvd_rho << rho
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/io/__init__.py", line 21, in __lshift__
self.write(u)
AttributeError: 'Conditional' object has no attribute '_cpp_object'
How do I correctly write density to pvd files? Thank you for your help!
