Function as_tensor expected in the definition of a scalar

variational formulation
1

Hi everyone,

I’m working with a phase field theory and I get when I try to derivate the local free energy function. It’s derivatives are sscalar values but I get an error saying Expecting a tuple of Index objects to A^indices := as_tensor(A, indices). The error points to the definition of f and the calling of the interpolation function h. Appreciate any help I can get.

Cheers

Here is my code:

from fenics import *

class InitialConditions(UserExpression):
    def eval(self, values, x):
        if near(x[0], lenght, DOLFIN_EPS) and near(x[1], height/2, DOLFIN_EPS): 
            values[0] = 0.0
        else:
            values[0] = 1.0
        values[1] = 1.0
    def value_shape(self):
        return (2,)

class CahnHilliardEquation(NonlinearProblem):
    def __init__(self, a, L, bcs):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
        self.bcs = bcs
    def F(self, b, x):
        assemble(self.L, tensor=b)
        for bc in self.bcs:
            bc.apply(b, x)
    def J(self, A, x):
        assemble(self.a, tensor=A)
        for bc in self.bcs:
            bc.apply(A)

class Gamma(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], lenght, DOLFIN_EPS) and near(x[1], height/2, DOLFIN_EPS)

# Output file
file = File("Results_AC/output.pvd", "compressed")

# Model parameters
A      = Constant(5.35e+07)
cse    = Constant(1.43e+02)
cle    = Constant(5.10e+00)
D      = Constant(8.50e-10)
M      = Constant( D / ( 2.0 *A ) )
L      = Constant(2.00e+00)
w      = Constant(1.03e+01)
alpha  = Constant(0.601281)

# Interpolation function
def h(n):
    return ( -2.0*n^3.0 + 3.0*n^2.0 )

# Double well potential
def g(n):
    return ( n^2.0 * ( 1.0 - n )^2.0 )

lmbda  = 1.0e02   # surface parameter
dt     = 5.0e-06  # time step
theta  = 1.0      # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson

# Geometry parameters and mesh parameters
lenght = 6.          # [mm]
height = 30.         # [mm]
Nelem  = 100         # [-]

# Form compiler options
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True

# Create mesh and build function space
mesh = RectangleMesh(Point(0.,  0.), Point(lenght, height), Nelem, Nelem*5, diagonal='right')
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1*P1)

# Define trial and test functions
du    = TrialFunction(ME)
q, v  = TestFunctions(ME)

# Define functions
u   = Function(ME)  # current solution
u0  = Function(ME)  # solution from previous converged step

# Split mixed functions
dc, dn = split(du)
c,  n  = split(u)
c0, n0 = split(u0)

# Create intial conditions and interpolate
u_init = InitialConditions(degree=1)
u.interpolate(u_init)
u0.interpolate(u_init)

# Compute the chemical potential df/dc
n = variable(n)
c = variable(c)
f    = A * ( c - cle - ( cse - cle ) * h(n) )^2.0 + w * g(n) 
dfdc = diff(f, c) 
dfdn = diff(f, n)

# Weak statement of the equations
L0 = ( c * q - c0 * q + dt * M * dot( grad( q ), grad( dfdc ) ) ) * dx
L1 = ( n * v - n0 * v + dt * L * alpha * dot( grad( v ), grad( n ) ) + dt * L * dfdn ) * dx
L = L0 + L1

# Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, u, du)

# Boundary conditions
bc1 = DirichletBC(ME.sub(0), Constant(0.0), Gamma(), method="pointwise")
bc2 = DirichletBC(ME.sub(1), Constant(0.0), Gamma(), method="pointwise")
bcs = [bc1, bc2]

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

# Step in time
t = 0.0
Nsteps = 50
T = Nsteps*dt
file << (u.split()[0], t)
while (t < T):
    t += dt
    u0.vector()[:] = u.vector()
    solver.solve(problem, u.vector())
    file << (u.split()[0], t)
2

Hi,
Several things strikes me when reading through your code.
First: The power operator is not ^, but **, you should change that throughout your code.
Second: The last part of L1 does not contain a test-function, and will lead to an arity-mismatch. Should probably be multiplied by v.

3

That was exactly it. Thanks for the help!