Implementing tanh() of fenics-function into variational formulation

variational formulation
1

Hello,

I’m having trouble implementing the following step function as a factor into my variational formulation:

\Theta (\mathbf{n}\cdot \mathbf{v}) = \frac{1}{2} (1 - \tanh(\frac{\mathbf{n}\cdot \mathbf{v}}{U_0 \delta}))

Where \mathbf{n} is the outward facing normal vector and \mathbf{v} is the 2D velocity vector. U_0 and \delta are constants. \mathbf{n} is always [1 \, 0]^T since the variational formulation is only applied at one specific vertical boundary.
\Theta should be then implemented into the variational form as follows:

F = \mathbf{n}\left[\frac{1}{2} \rho |\mathbf{v}|^2 \Theta (\mathbf{n}\cdot \mathbf{v}) \right]

Since I’m multiplying this term with ds(1) in FEniCS, I’m making sure that this is only applied at the specific boundary marked with 1.

The whole boundary term is only relevant when the velocity x-component is negative, then \Theta becomes 0.5, otherwise it’s 0. \Theta is basically a conditional statement.

In my code \mathbf{v} is a function of a 2D function space: v = Function(FunctionSpace) And therefore a list tensor. \mathbf{n} is defined as n = FacetNormal(mesh).

This minimal code snippet reproduces my error:

from fenics import *
from ufl import tanh

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0.0)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 96)

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0.0)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 96)

def boundaries(x, on_boundary):
    return on_boundary

# Define Function Space
mesh = UnitSquareMesh(96, 96)

P2v = VectorElement("Lagrange", mesh.ufl_cell(), 2)

ME = FunctionSpace(mesh, P2v)

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)

# mark subdomains to boundaries
Right().mark(boundaries, 1)
Left().mark(boundaries, 2)
Top().mark(boundaries,3)
Bottom().mark(boundaries,4)

# Define boundary measure on Neumann part of boundary
ds = Measure("ds", subdomain_data=boundaries)
n = FacetNormal(mesh)

bcv_top                 = DirichletBC(ME, Constant((0,0)) , boundaries, 3)
bcv_bottom              = DirichletBC(ME, Constant((0,0)) , boundaries, 4)
bcv_inflow              = DirichletBC(ME, Constant((1,0)), boundaries, 2)


# put all velocity boundaries in list
bcs = [bcv_inflow, bcv_top, bcv_bottom]

# Define trial and test functions

dv       = TrialFunction(ME)
w        = TestFunction(ME)

# # Define functions
v   = Function(ME)
v0  = Function(ME)

rho     = Constant(924.3)
U       = 0.0048
delta   = 0.0025

theta = 0.5 * (1 - tanh(dot(v,n)/(U*delta)))
F = 0.5 * rho * abs(v0)**2 * theta * dot(n, w) * ds(1)
a, L = lhs(F), rhs(F)


problem = LinearVariationalProblem(a, L, v, bcs=bcs)
solver  = LinearVariationalSolver(problem)

My error message is:

Form has no parts with arity 2.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Expecting scalar.
Traceback (most recent call last):
  File "/Users/pauline/Desktop/minimal code.py", line 68, in <module>
    problem = LinearVariationalProblem(a, L, v, bcs=bcs)
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/dolfin/fem/problem.py", line 55, in __init__
    L = Form(L, form_compiler_parameters=form_compiler_parameters)
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/dolfin/fem/form.py", line 43, in __init__
    ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/dolfin/jit/jit.py", line 47, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/dolfin/jit/jit.py", line 97, in ffc_jit
    return ffc.jit(ufl_form, parameters=p)
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/ffc/jitcompiler.py", line 217, in jit
    module = jit_build(ufl_object, module_name, parameters)
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/ffc/jitcompiler.py", line 130, in jit_build
    module, signature = dijitso.jit(jitable=ufl_object,
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/dijitso/jit.py", line 165, in jit
    header, source, dependencies = generate(jitable, name, signature, params["generator"])
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/ffc/jitcompiler.py", line 65, in jit_generate
    code_h, code_c, dependent_ufl_objects = compile_object(ufl_object,
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/ffc/compiler.py", line 142, in compile_form
    return compile_ufl_objects(forms, "form", object_names,
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/ffc/compiler.py", line 190, in compile_ufl_objects
    ir = compute_ir(analysis, prefix, parameters, jit)
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/ffc/representation.py", line 182, in compute_ir
    irs = [_compute_integral_ir(fd, form_id, prefix, element_numbers, classnames, parameters, jit)
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/ffc/representation.py", line 182, in <listcomp>
    irs = [_compute_integral_ir(fd, form_id, prefix, element_numbers, classnames, parameters, jit)
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/ffc/representation.py", line 455, in _compute_integral_ir
    ir = r.compute_integral_ir(itg_data,
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/ffc/uflacs/uflacsrepresentation.py", line 132, in compute_integral_ir
    uflacs_ir = build_uflacs_ir(itg_data.domain.ufl_cell(),
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/ffc/uflacs/build_uflacs_ir.py", line 338, in build_uflacs_ir
    V, V_deps, V_targets = build_scalar_graph(expressions)
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/ffc/uflacs/build_uflacs_ir.py", line 823, in build_scalar_graph
    scalar_expressions = rebuild_with_scalar_subexpressions(G)
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/ffc/uflacs/analysis/graph_rebuild.py", line 276, in rebuild_with_scalar_subexpressions
    ws = reconstruct_scalar_subexpressions(v, wops)
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/ufl/corealg/multifunction.py", line 100, in __call__
    return self._handlers[o._ufl_typecode_](o, *args)
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/ffc/uflacs/analysis/graph_rebuild.py", line 62, in scalar_nary
    error("Expecting scalar.")
  File "<string>", line 1, in <lambda>
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/ufl/log.py", line 172, in error
    raise self._exception_type(self._format_raw(*message))
Exception: Expecting scalar.

I have the “Form has no parts of arity 2” error because the variational form isn’t complete. I don’t have this error in the complete code. But what’s important is the “expecting scalar” error message.

How can I implement this step function? Does anybody have any idea?
Kind regards,
Pauline

2

Hi,
when running your code I got the error that ‘states’ in

problem = LinearVariationalProblem(a, L, states, bcs=bcs)

is not defined.

3

F is vector-valued, because abs(v0) is the absolute value of the (vector-valued) function v0, not the norm/magnitude of v0. Assuming the notation

\left| \mathbf{v} \right|^2

refers to the squared magnitude of \mathbf{v}, i.e.

\left| \mathbf{v} \right|^2 = \mathbf{v} \cdot \mathbf{v}

you could implement this as:

F = 0.5 * rho * dot(v0,v0) * theta * dot(n, w) * ds(1)
1 Like
4

Hi, thanks for the quick reply. I fixed the typo.

1 Like
5

Thanks for the quick reply. This fixed my problem!