ufl.log.UFLException: Product can only represent products of scalar

Question

Hi, I have added a part (as per the advection-diffusion demo) to my code to alleviate numerical instability. Yet I get this error message which I can't figure out where it relates to. Could you please help? Many thanks,

(parts of the error message have been deleted to fit the space)

Product can only represent products of scalars.
FAILURE in reuse_if_possible: type(o) = <class 'ufl.algebra.Product'> operands =

w_{10}

-0.0864 + -0.1728 + 0.864 * w_2 / (0.0025 + w_2)



<class 'ufl.algebra.Product'> ; (sum_{i_{20}} ((w_5)[i_{20}]) * (({ A

<class 'ufl.algebra.Sum'> ; 0.02 * (-1 * (v_{-1} * (-0.0864 + -0.1728


Traceback (most recent call last):   File "s13.py", line 240, in <module>
    F = action(F, S)   File "C:\FEniCS\lib\site-packages\ufl\formoperators.py", line 100,
    return compute_form_action(form, coefficient)   File "C:\FEniCS\lib\site-packages\ufl\algorithms\formtransformation  389, in compute_form_action
    return replace(form, { u: coefficient })   File "C:\FEniCS\lib\site-packages\ufl\algorithms\transformations.py , in replace
    return apply_transformer(e, Replacer(mapping2))   File "C:\FEniCS\lib\site-packages\ufl\algorithms\transformations.py , in apply_transformer
    return transform_integrands(e, _transform)   File "C:\FEniCS\lib\site-packages\ufl\algorithms\transformations.py , in transform_integrands
    integrand = transform(itg.integrand())   File "C:\FEniCS\lib\site-packages\ufl\algorithms\transformations.py , in
_transform
    return transformer.visit(expr)   File "C:\FEniCS\lib\site-packages\ufl\algorithms\transformations.py , in visit
    r = h(o, *map(self.visit, o.operands()))   File "C:\FEniCS\lib\site-packages\ufl\algorithms\transformations.py , in visit
    r = h(o, *map(self.visit, o.operands()))   File "C:\FEniCS\lib\site-packages\ufl\algorithms\transformations.py , in visit
    r = h(o, *map(self.visit, o.operands()))   File "C:\FEniCS\lib\site-packages\ufl\algorithms\transformations.py , in visit
    r = h(o, *map(self.visit, o.operands()))   File "C:\FEniCS\lib\site-packages\ufl\algorithms\transformations.py , in visit
    r = h(o, *map(self.visit, o.operands()))   File "C:\FEniCS\lib\site-packages\ufl\algorithms\transformations.py , in visit
    r = h(o, *map(self.visit, o.operands()))   File "C:\FEniCS\lib\site-packages\ufl\algorithms\transformations.py , in visit
    r = h(o, *map(self.visit, o.operands()))   File "C:\FEniCS\lib\site-packages\ufl\algorithms\transformations.py , in visit
    r = h(o, *map(self.visit, o.operands()))   File "C:\FEniCS\lib\site-packages\ufl\algorithms\transformations.py , in reuse_if_possible
    r = o.reconstruct(*operands)   File "C:\FEniCS\lib\site-packages\ufl\expr.py", line 214, in recons
    return self.__class__._uflclass(*operands)   File "C:\FEniCS\lib\site-packages\ufl\algebra.py", line 177, in __n
    error("Product can only represent products of scalars.")   File "C:\FEniCS\lib\site-packages\ufl\log.py", line 148, in error
    raise self._exception_type(self._format_raw(*message)) ufl.log.UFLException: Product can only represent products of scalars.

Comments

The code that faces this error message is the following (beginning parts deleted to fit space):

The source function f is nonlinear in the variable, hence calling the nonlinear solver.
S and P are variables determined from separate PDEs earlier in the code and fS is a defined function of S:

fS = (B*S/(K+S)) - C 

(B, C, K are constants)

def boundary_value(n):
    if n < 10:
        return float(n)/10.0
    else:
        return 1.0

mesh = Rectangle(0,0,1,H,divx,divy)
subdomains = MeshFunction("uint", mesh, 2)
class Biomass(SubDomain):
    def inside(self, x, on_boundary):
            return True if x[1] <= 0.2 + 0.05*sin(4*pi*x[0]) else False
class Interface(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[1] >= 0.2 + 0.05*sin(4*pi*x[0]) else False

subdomain0 = Biomass()
subdomain0.mark(subdomains, 0)
subdomain1 = Interface()
subdomain1.mark(subdomains, 1)

h = CellSize(mesh)

# Create FunctionSpaces
Q = FunctionSpace(mesh, "CG", 1)
V = VectorFunctionSpace(mesh, "CG", 2)
vel = -grad(P)
velocity = project(vel, V)

u0 = Function(Q)
# Parameters
T = 0.1
dt = 0.02
t = dt
c = 0.00005

# Test and trial functions
u, v = TrialFunction(Q), TestFunction(Q)

f = u*(fS - Dl) - u*u*fS

# Mid-point solution
u_mid = 0.5*(u0 + u)

# Residual
r = u-u0 + dt*(dot(velocity, grad(u_mid)) - c*div(grad(u_mid)) - f)

# Galerkin variational problem
F = v*(u-u0)*dx + dt*(v*dot(velocity, grad(u_mid))*dx + c*dot(grad(v), grad(u_mid))*dx)

# Add SUPG stabilisation terms
vnorm = sqrt(dot(velocity, velocity))
F += (h/2.0*vnorm)*dot(velocity, grad(v))*r*dx

# Set up boundary condition
g = Constant(1.0)
bcs = DirichletBC(Q, g, Interface())
S = Function(V)  
out_file = File("results.pvd")
# Set intial condition
u = u0
# time-stepping
while t < T:
    F = action(F, S)
    J = derivative(F, S, u)
    # compute solution
    problem = NonlinearVariationalProblem(F, S, bcs, J)
    solver  = NonlinearVariationalSolver(problem)
    solver.solve()
    u0 = u
    plot(u)
    out_file << (u, t)
    t += dt
    g.assign( boundary_value( int(t/dt) ) )
interactive()

---

You need to make your code much simpler but complete (i.e., runnable) to get an answer. You can't expect others to debug your whole solver.

---

Sorry. That certainly wasn't the intention.
I shall definitely abide by this next time.

Answer 1

F = action(F, S)

that's nonsense. F is (non-linear) map Q$\times$Q$\rightarrow{\cal R}$ (which is linear in second argument) and you want to apply it on function from V.

Comments

Understood. My bad.