Exception: Expecting scalar

variational formulation
1

Hi,

Using dolfin-2019.1. Here is the relevant code excerpt:

import dolfin as dlf
DT = dlf.Constant(1e-4)  # time step
MU1 = dlf.Constant(1.0e-3)  # dummy dynamic viscosity value
RHO1 = dlf.Constant(917)  # density
G = dlf.Constant((0, -9.81))

# 2D rectangular channel mesh ----------
mesh = dlf.RectangleMesh(dlf.Point(0, 0), dlf.Point(0.5, 2.0), 10, 40, 'right')
n = dlf.FacetNormal(mesh)

# function spaces ----------------------
V = dlf.VectorFunctionSpace(mesh, 'DG', 1)
Q = dlf.FunctionSpace(mesh, 'CG', 2)
# trial and test functions
u = dlf.TrialFunction(V)  # phase velocity trial function
v = dlf.TestFunction(V)
jump_v = v('+') - v('-')

p = dlf.TrialFunction(Q)  # shared pressure field for both phases
vf = dlf.TrialFunction(Q)  # phase volume fraction trial function
q = dlf.TestFunction(Q)

u1_o = dlf.Function(V); u1_o.assign(dlf.Constant((0, 1.0)))
u1_tld = dlf.Function(V); u1_tld.assign(dlf.Constant((0, 1.0)))
u1_str = dlf.Function(V); u1_str.assign(dlf.Constant((0, 1.0)))

u_r_tld = u1_tld - u0_tld

# pressures
p_o = dlf.Function(Q); p_o.assign(dlf.Constant(0.0))
p_tnt = dlf.Function(Q); p_tnt.assign(dlf.Constant(0.0))
p_str = dlf.Function(Q); p_str.assign(dlf.Constant(0.0))
delta_p = dlf.Function(Q); delta_p.assign(dlf.Constant(0.0))

# volume fractions
vf0_o = dlf.Function(Q); vf0_o.assign(dlf.Constant(0.95))
vf0_tld = dlf.Function(Q); vf0_tld.assign(dlf.Constant(0.95))
vf0_tnt = dlf.Function(Q); vf0_tnt.assign(dlf.Constant(0.95))
#
vf1_o = dlf.Function(Q); vf1_o.assign(dlf.Constant(0.05))
vf1_tld = dlf.Function(Q); vf1_tld.assign(dlf.Constant(0.05))
vf1_tnt = dlf.Function(Q); vf1_tnt.assign(dlf.Constant(0.05))

# continuity equation of phase 1
#
N1 = (1/DT) * dlf.inner(q, vf - vf1_o) * dlf.dx
#
c_C1 = vf1_tld * dlf.inner(u1_tld, dlf.grad(q)) * dlf.dx
#
# volume fractions surface integrals,
# OMITTING THE - LIKE C0 & C1, TO ACCOUNT FOR IT IN THE EQNS
#
r1 = dlf.inner(vf1_tld, q) * dlf.dot(u1_tld, n) * dlf.ds
r0 = dlf.inner(vf0_tld, q) * dlf.dot(u0_tld, n) * dlf.ds
#
cont1 = N1 - c_C1 + r1
c_lhs1 = dlf.lhs(cont1)
c_rhs1 = dlf.rhs(cont1)
dlf.solve(c_lhs1 == c_rhs1, vf1_tnt, dbcvf1)
vf1_tld = THETA_NL * vf1_tnt + (1 - THETA_NL) * vf1_o

# velocity prediction
#
# discrete phase (1)
u_n_1 = 0.5 * (dlf.dot(u1_tld, n) + abs(dlf.dot(u1_tld, n)))
gradu_n_1 = 0.5 * (dlf.dot(dlf.grad(u1_tld), n)
                   + abs(dlf.dot(dlf.grad(u1_tld), n)))
#
M1 = (RHO1 * vf1_tld / DT) * dlf.inner(v, u - u1_o) * dlf.dx
#
jump_A1 = dlf.inner(u_n_1('+') * u('+') - u_n_1('-') * u('-'), jump_v) \
        * vf1_tld * RHO1 * dlf.dS
A1 = - RHO1 * vf1_tld * dlf.inner(dlf.grad(v) * u1_tld, u) * dlf.dx \
     - dlf.div(vf1_tld * RHO1 * u1_tld) * dlf.inner(v, u) * dlf.dx \
     + vf1_tld * RHO1 * dlf.dot(n, u1_tld) * dlf.inner(v, u) * dlf.ds
#
jump_K1 = dlf.inner(gradu_n_1('+') - gradu_n_1('-'), jump_v) \
        * vf1_tld * MU1 * dlf.dS
K1 = vf1_tld * MU1 * dlf.inner(dlf.grad(v), dlf.grad(u)) * dlf.dx \
     - vf1_tld * MU1 * dlf.inner(dlf.dot(n, dlf.grad(u1_tld)), v) * dlf.ds
#
F1 = 0.75 * C * vf1_tld * vf0_tld * RHO0 * dlf.inner(v, u) \
     * dlf.sqrt (dlf.dot(u_r_tld, u_r_tld)) * dlf.dx
#
# C1 = - C1 from paper to avoid bad operand type for unary error
# ADD THE - BACK IN THE pred1 FORM BELOW
C1 = vf1_tld * dlf.inner(v, dlf.grad(p_str)) * dlf.dx
#
b1 = vf1_tld * RHO1 *dlf.inner(G, v) * dlf.dx
#
f1 = 0.75 * C * vf1_tld * vf0_tld * RHO0 * dlf.inner(v, u0_tld) \
     * dlf.sqrt(dlf.dot(u_r_tld, u_r_tld)) * dlf.dx
#
pred1 = M1 + A1 + jump_A1 + K1 + jump_K1 + F1 + C1 - b1 - f1
lhs1 = dlf.lhs(pred1)
rhs1 = dlf.rhs(pred1)
dlf.solve(lhs1 == rhs1, u1_str, dbcu)

At the last line I get the error below where a scalar is expected somewhere in pred1 and I cannot figure out where it comes from. Anyone can help?

python3 two_phase_FEM.py 
created mesh with  800 cells and  451 vertices
Solving linear variational problem.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Expecting scalar.
Traceback (most recent call last):
  File "two_phase_FEM.py", line 165, in <module>
    dlf.solve(lhs1 == rhs1, u1_str, dbcu)
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/dolfin/fem/solving.py", line 220, in solve
    _solve_varproblem(*args, **kwargs)
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/dolfin/fem/solving.py", line 242, in _solve_varproblem
    form_compiler_parameters=form_compiler_parameters)
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/dolfin/fem/problem.py", line 55, in __init__
    L = Form(L, form_compiler_parameters=form_compiler_parameters)
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/dolfin/fem/form.py", line 44, in __init__
    mpi_comm=mesh.mpi_comm())
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/dolfin/jit/jit.py", line 47, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/dolfin/jit/jit.py", line 97, in ffc_jit
    return ffc.jit(ufl_form, parameters=p)
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/ffc/jitcompiler.py", line 217, in jit
    module = jit_build(ufl_object, module_name, parameters)
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/ffc/jitcompiler.py", line 133, in jit_build
    generate=jit_generate)
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/dijitso/jit.py", line 165, in jit
    header, source, dependencies = generate(jitable, name, signature, params["generator"])
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/ffc/jitcompiler.py", line 66, in jit_generate
    prefix=module_name, parameters=parameters, jit=True)
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/ffc/compiler.py", line 143, in compile_form
    prefix, parameters, jit)
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/ffc/compiler.py", line 190, in compile_ufl_objects
    ir = compute_ir(analysis, prefix, parameters, jit)
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/ffc/representation.py", line 183, in compute_ir
    for (form_id, fd) in enumerate(form_datas)]
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/ffc/representation.py", line 183, in <listcomp>
    for (form_id, fd) in enumerate(form_datas)]
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/ffc/representation.py", line 460, in _compute_integral_ir
    parameters)
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/ffc/uflacs/uflacsrepresentation.py", line 139, in compute_integral_ir
    parameters)
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/ffc/uflacs/build_uflacs_ir.py", line 338, in build_uflacs_ir
    V, V_deps, V_targets = build_scalar_graph(expressions)
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/ffc/uflacs/build_uflacs_ir.py", line 823, in build_scalar_graph
    scalar_expressions = rebuild_with_scalar_subexpressions(G)
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/ffc/uflacs/analysis/graph_rebuild.py", line 276, in rebuild_with_scalar_subexpressions
    ws = reconstruct_scalar_subexpressions(v, wops)
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/ufl/corealg/multifunction.py", line 100, in __call__
    return self._handlers[o._ufl_typecode_](o, *args)
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/ffc/uflacs/analysis/graph_rebuild.py", line 62, in scalar_nary
    error("Expecting scalar.")
  File "<string>", line 1, in <lambda>
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/ufl/log.py", line 172, in error
    raise self._exception_type(self._format_raw(*message))
Exception: Expecting scalar.
2

Can you turn this into a minimal working example? There are many missing variable definitions.

3

Here’s a mwe up to the point where the error occurs. I tested it and it runs identically to the full code with the same error.

#!/usr/bin/python3
""" Two-phase flow transport model based on Gosman model with P1DG-P2
    element (Cotter et al. Journal of Computational Physics 228, 2009),
    and the numerical approach of Jacobs et al. (Geophysical Journal
    International 192 (2013)).
"""

# import time
# import numpy as np
import dolfin as dlf
# import meshio as mio

# numerical parameters -----------------
DT = dlf.Constant(1e-4)  # time step
NUM_STEPS = 100  # number of time steps
N_SAVE = 20  # solution saving interval
T_FINAL = NUM_STEPS * DT  # final solution time
THETA_NL = dlf.Constant(0.5)  # Picard iteration relaxation parameter

# phase properties ---------------------
# continuous phase 0 (air)
MU0 = dlf.Constant(1.63e-5)  # dynamic viscosity
RHO0 = dlf.Constant(1.38)  # density
NU0 = MU0 / RHO0  # kinematic viscosity

# discrete phase 1 (snow)
MU1 = dlf.Constant(1.0e-3)  # dummy dynamic viscosity value
RHO1 = dlf.Constant(917)  # density
NU1 = MU1 / RHO1  # kinematic viscosity
D1 = dlf.Constant(1e-4)  # particle diameter
C = dlf.Constant(0.2)  # drag coefficient REWRITE AS FUNCTION OF PARTICLE REYNOLDS
G = dlf.Constant((0, -9.81))

# 2D rectangular channel mesh ----------
mesh = dlf.RectangleMesh(dlf.Point(0, 0), dlf.Point(0.5, 2.0),
                         10, 40, 'right')
NUM_CELLS = mesh.num_cells()
NUM_VERTICES = mesh.num_vertices()
print("created mesh with {:4d} cells and {:4d} vertices"
      .format(NUM_CELLS, NUM_VERTICES))
dlf.File('mesh.pvd').write(mesh)
n = dlf.FacetNormal(mesh)

# Boundaries ---------------------------
# need x = dlf.SpatialCoordinate(mesh)? Doesn't look like it anymore.
LEFT = 'near(x[1], 0)'
RIGHT = 'near(x[1], 0.5)'
INLET = 'near(x[0], 0)'
OUTLET = 'near(x[1], 2.0)'

# function spaces ----------------------
V = dlf.VectorFunctionSpace(mesh, 'DG', 1)
Q = dlf.FunctionSpace(mesh, 'CG', 2)

# Boundary conditions ------------------
# velocities
dbcu_left = dlf.DirichletBC(V, dlf.Constant((0.0, 0.0)), LEFT,
                           method='geometric')
dbcu_right = dlf.DirichletBC(V, dlf.Constant((0.0, 0.0)), RIGHT,
                              method='geometric')
dbcu_inlet = dlf.DirichletBC(V, dlf.Constant((0.0, 1.0)), INLET,
                             method='geometric')
dbcu = [dbcu_left, dbcu_right, dbcu_inlet]

# pressure
dbcp_outlet = dlf.DirichletBC(Q, dlf.Constant(0.0), OUTLET)
dbcp = [dbcp_outlet]

# discrete phase volume fraction
dbcvf1_inlet = dlf.DirichletBC(Q, dlf.Constant(0.05), INLET)
dbcvf1 = [dbcvf1_inlet]

# trial and test functions
u = dlf.TrialFunction(V)  # phase velocity trial function
v = dlf.TestFunction(V)
jump_v = v('+') - v('-')

p = dlf.TrialFunction(Q)  # shared pressure field for both phases
vf = dlf.TrialFunction(Q)  # phase volume fraction trial function
q = dlf.TestFunction(Q)

# old time (*_o), tilde (*_tld) and current time (*_str) fields
# velocities
u0_o = dlf.Function(V); u0_o.assign(dlf.Constant((0, 1.0)))
u0_tld = dlf.Function(V); u0_tld.assign(dlf.Constant((0, 1.0)))
u0_str = dlf.Function(V); u0_str.assign(dlf.Constant((0, 1.0)))
#
u1_o = dlf.Function(V); u1_o.assign(dlf.Constant((0, 1.0)))
u1_tld = dlf.Function(V); u1_tld.assign(dlf.Constant((0, 1.0)))
u1_str = dlf.Function(V); u1_str.assign(dlf.Constant((0, 1.0)))
#
u_r_tld = u1_tld - u0_tld

# pressures
p_o = dlf.Function(Q); p_o.assign(dlf.Constant(0.0))
p_tnt = dlf.Function(Q); p_tnt.assign(dlf.Constant(0.0))
p_str = dlf.Function(Q); p_str.assign(dlf.Constant(0.0))
delta_p = dlf.Function(Q); delta_p.assign(dlf.Constant(0.0))

# volume fractions
vf0_o = dlf.Function(Q); vf0_o.assign(dlf.Constant(0.95))
vf0_tld = dlf.Function(Q); vf0_tld.assign(dlf.Constant(0.95))
vf0_tnt = dlf.Function(Q); vf0_tnt.assign(dlf.Constant(0.95))
#
vf1_o = dlf.Function(Q); vf1_o.assign(dlf.Constant(0.05))
vf1_tld = dlf.Function(Q); vf1_tld.assign(dlf.Constant(0.05))
vf1_tnt = dlf.Function(Q); vf1_tnt.assign(dlf.Constant(0.05))


# velocity prediction
#
# discrete phase (1)
u_n_1 = 0.5 * (dlf.dot(u1_tld, n) + abs(dlf.dot(u1_tld, n)))
gradu_n_1 = 0.5 * (dlf.dot(dlf.grad(u1_tld), n)
                   + abs(dlf.dot(dlf.grad(u1_tld), n)))
#
M1 = (RHO1 * vf1_tld / DT) * dlf.inner(v, u - u1_o) * dlf.dx
#
jump_A1 = dlf.inner(u_n_1('+') * u('+') - u_n_1('-') * u('-'), jump_v) \
        * vf1_tld * RHO1 * dlf.dS
A1 = - RHO1 * vf1_tld * dlf.inner(dlf.grad(v) * u1_tld, u) * dlf.dx \
     - dlf.div(vf1_tld * RHO1 * u1_tld) * dlf.inner(v, u) * dlf.dx \
     + vf1_tld * RHO1 * dlf.dot(n, u1_tld) * dlf.inner(v, u) * dlf.ds
#
jump_K1 = dlf.inner(gradu_n_1('+') - gradu_n_1('-'), jump_v) \
        * vf1_tld * MU1 * dlf.dS
K1 = vf1_tld * MU1 * dlf.inner(dlf.grad(v), dlf.grad(u)) * dlf.dx \
     - vf1_tld * MU1 * dlf.inner(dlf.dot(n, dlf.grad(u1_tld)), v) * dlf.ds
#
F1 = 0.75 * C * vf1_tld * vf0_tld * RHO0 * dlf.inner(v, u) \
     * dlf.sqrt (dlf.dot(u_r_tld, u_r_tld)) * dlf.dx
#
# C1 = - C1 from paper to avoid bad operand type for unary error
# ADD THE - BACK IN THE pred1 FORM BELOW
C1 = vf1_tld * dlf.inner(v, dlf.grad(p_str)) * dlf.dx
#
b1 = vf1_tld * RHO1 *dlf.inner(G, v) * dlf.dx
#
f1 = 0.75 * C * vf1_tld * vf0_tld * RHO0 * dlf.inner(v, u0_tld) \
     * dlf.sqrt(dlf.dot(u_r_tld, u_r_tld)) * dlf.dx
#
pred1 = M1 + A1 + jump_A1 + K1 + jump_K1 + F1 + C1 - b1 - f1
lhs1 = dlf.lhs(pred1)
rhs1 = dlf.rhs(pred1)
dlf.solve(lhs1 == rhs1, u1_str, dbcu)
4

This term doesn’t make sense:

abs(dlf.dot(u1_tld, n))

I’m not familiar with the method you’re wanting to use, so depending on what you formulation needs consider taking the absolute value of each component

dlf.as_vector(list(map(abs, dlf.dot(dlf.grad(u1_tld), n))))
5

Using your component by component suggestion results in the following shape mismatch error

Can't add expressions with different shapes.
Traceback (most recent call last):
  File "mwe.py", line 114, in <module>
    u_n_1 = 0.5 * (dlf.dot(u1_tld, n) + dlf.as_vector(list(map(abs, dlf.dot(dlf.grad(u1_tld), n)))))
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/ufl/exproperators.py", line 212, in _add
    return Sum(self, o)
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/ufl/algebra.py", line 53, in __new__
    error("Can't add expressions with different shapes.")
  File "/home/zbinkz/HGST/Projects/Python/FEniCS-p3/lib/python3.7/site-packages/ufl/log.py", line 172, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Can't add expressions with different shapes.

The upwinding expression I use worked for a scalar DG advection formulation I used for another code with FEniCS and is not reporting errors here if I comment the solve statement at the end. The scalar error here also appears with LinearVariational* classes. Anyway but I will try another jump function to see if it goes away.

BTW I am trying to reproduce the Euler-Euler approach described here:
Geophys. J. Int. (2013) 192, 647–665, doi: 10.1093/gji/ggs059
Multiphase flow modelling of volcanic ash particle settling in water using adaptive unstructured meshes
C. T. Jacobs, G. S. Collins, M. D. Piggott, S. C. Kramer and C. R. G. Wilson

Would you recommend dolfin_dg?

6

abs expects a scalar since it computes the absolute value. Perhaps you’re expecting it to give you a norm of the vector, something like ufl.sqrt(ufl.dot(u, u))?

I’m not familiar with the paper you’ve posted, but well aware of fluidity’s popularity for such problems. I can’t really recommend the best tool for what you want to achieve since you know your requirements much better than I do.

7

abs is not the issue here. There is something else in the pred1 form that is causing the error and I guess I will need to find it through term by term elimination.

I’m fairly new to dolfin, ufl and DG so that is why I asked about dolfin_dg. Would it make life easier to a novice trying to code algorithms of the type in Fluidity? Of course I don’t expect you to be familiar with the inner workings of Fluidity but just asking.

8

I found the culprit:

jump_K1 = dlf.inner(gradu_n_1('+') - gradu_n_1('-'), jump_v) * vf1_tld * MU1 * dlf.dS

This is my botched attempt at a jump function for the inner element surface integral of the stiffness matrix term K1 surface integral resulting from the IBP:

vf1_tld * MU1 * dlf.inner(dlf.dot(n, dlf.grad(u1_tld)), v) * dlf.ds

Here, vf1_tld is a scalar-valued continuous function, MU1 is a Constant(), n is the face normal, u1_tld is a vector-valued discontinuous function and v is the vector-valued discontinuous test function. So the real question becomes how to define a half-decent jump function for K1. Any ideas? It somewhere in the Fluidity code but a tip would be highly appreciated :grin: