cahn-hilliard-navier-stokes:ValueError: not enough values to unpack (expected 4, got 2)

variational formulation
1

I am trying to couple NS equation and CH equation, I have some problems, here are my code and error message, I hope you can help me solve this problem.

import random
from numpy.ma import tanh
from fenics import *  



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
        values[2] = (0.0, 0.0)
        values[3] = 0.0

    def value_shape(self):
        return 4,


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


def rho(c_inp):
    return conditional(lt(c_inp, c_ox), left(c_inp), conditional(lt(c_inp, c_met), mid(c_inp), right(c_inp)))



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
yita = 10
g = (0.0, 9.80665)

mesh = RectangleMesh(Point(0, 0), Point(L, 2 * L), nx, ny)

P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
ME = FunctionSpace(mesh, P1 * P1 * P2 * P1)  # c mu [v0,v1],p


d_c_mu_u_p = TrialFunction(ME)
fai, pot, v, q = TestFunctions(ME)
c_mu_u_p = Function(ME)  # current solution
c_mu_u_p_0 = Function(ME)  # solution from previous converged step

dc, dmu, du, dp = split(d_c_mu_u_p)
c, mu, u, p = split(c_mu_u_p)
c0, mu0, u0, p0 = split(c_mu_u_p_0)


c_mu_init = InitialConditions(degree=1)
c_mu_u_p.interpolate(c_mu_init)
c_mu_u_p_0.interpolate(c_mu_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 * fai * dx - c0 * fai * dx + dt * dot(dot(u, grad(c)), fai) * dx + dt * M * dot(grad(mu), grad(fai)) * dx
L1 = mu * pot * dx - df_dc * pot * dx - coef * dot(grad(c), grad(pot)) * dx
L2 = dot(u, v) * dx - \
     dot(u0, v) * dx + \
     dt * dot(grad(u) * u, v) * dx + \
     dt / rho_int_ox * dot(grad(p), v) * dx + \
     yita * dt / rho_int_ox * dot(grad(u), grad(v)) * dx + \
     dt / rho_int_ox * dot(c * grad(mu), v) * dx + \
     dt / rho_int_ox * dot(rho(c) * g, v) * dx
L3 = div(u) * q*dx
L = L0 + L1 + L2 + L3

a = derivative(L, c_mu_u_p, d_c_mu_u_p)

problem = CahnHilliardEquation(a, L)
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6

t = 0.0
solver.solve(problem, c_mu_u_p.vector())

Here is the error message:

Traceback (most recent call last):
  File "/home/wangyao/PycharmProjects/hello_word/main.py", line 114, in <module>
    fai, pot, v, q = TestFunctions(ME)
ValueError: not enough values to unpack (expected 4, got 2)

Thank you very much!

The governing equation is as follows:
image

I believe the following is the simplest code:

from fenics import *

L = 0.176
nx = 192
ny = 384

mesh = RectangleMesh(Point(0, 0), Point(L, 2 * L), nx, ny)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
ME = FunctionSpace(mesh, P1 * P1 * P2 * P1)  # c mu [v0,v1],p

d_c_mu_u_p = TrialFunction(ME)
fai, pot, v, q = TestFunctions(ME)
c_mu_u_p = Function(ME)  # current solution
c_mu_u_p_0 = Function(ME)  # solution from previous converged step

dc, dmu, du, dp = split(d_c_mu_u_p)
c, mu, u, p = split(c_mu_u_p)
c0, mu0, u0, p0 = split(c_mu_u_p_0)
2

This is because you are using the product operator between the elements, instead of usinga mixed element:

from fenics import *  
L = 1
nx = 5
ny = 5
mesh = RectangleMesh(Point(0, 0), Point(L, 2 * L), nx, ny)

P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
element = MixedElement([P1, P1, P2, P1])
ME = FunctionSpace(mesh, element)

d_c_mu_u_p = TrialFunction(ME)
fai, pot, v, q = TestFunctions(ME)

Next time, please read the guidelines: Read before posting: How do I get my question answered? - #3 before posting a question, as it is quite clear that your error has nothing to do with your PDE, and that you could reduce your example to ~10 lines, as I’ve done above.

3

Thank you dokken, there is still some confusion because I noticed that the element product operation is used in the fenapack NS equation example,

# Build Taylor-Hood function space
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, P2*P1)

So why doesn’t the product of this finite element apply to my case?
link:https://github.com/blechta/fenapack/blob/master/demo/navier-stokes-pcd/demo_navier-stokes-pcd.py

4

Because the product is between only two elements.
When using products, you create a mixed-element per product, which can be seen by printing:

print( P1*P1*P2*P1)

which results in

<Mixed element: (<Mixed element: (<Mixed element: (<CG1 on a triangle>, <CG1 on a triangle>)>, <vector element with 2 components of <CG2 on a triangle>>)>, <CG1 on a triangle>)>

while my proposed syntax produces a single mixed element

print( MixedElement([P1, P1, P2, P1]))                                  

<Mixed element: (<CG1 on a triangle>, <CG1 on a triangle>, <vector element with 2 components of <CG2 on a triangle>>, <CG1 on a triangle>)>
1 Like
5

Thank you dokken and it’s very useful.