Issue with Error at the Program

Errors
1

I am encountering an error while executing this program. The issue arises specifically at the last line. Could you please provide feedback or suggestions to resolve it?

from fenics import *
import numpy as np

channel_length, channel_height = 1, 0.5
mesh = RectangleMesh(Point(0, 0), Point(channel_length, channel_height), 100, 100) 

h = Constant(mesh.hmax())
n = FacetNormal(mesh)

Q = FunctionSpace(mesh, "DG", 1)

p_trial = TrialFunction(Q)
q = TestFunction(Q)
p_prev = Function(Q)
p = Function(Q)

u = Constant((0,0))

class InletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - 0) < DOLFIN_EPS and on_boundary

class OutletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - channel_length) < DOLFIN_EPS and on_boundary

inlet = InletBoundary()
outlet = OutletBoundary()

class DirichletBoundary_other(SubDomain):
    def inside(self, x, on_boundary):
        return (
            on_boundary
            and not (
                inlet.inside(x, on_boundary) or
                outlet.inside(x, on_boundary)
            )
        )

other = DirichletBoundary_other()

dx = Measure('dx', domain=mesh)
ds = Measure('ds', domain=mesh)                                                 
dS = Measure('dS', domain=mesh)                                                 

boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_markers.set_all(0)
other.mark(boundary_markers, 1)
inlet.mark(boundary_markers, 2)
outlet.mark(boundary_markers, 3)

ds = Measure("ds", domain=mesh, subdomain_data=boundary_markers)

dt = 0.001
beta = 500
beta_2 = 5000

b_volume = -inner(grad(p_prev - p_trial), grad(q)) * dx + 1/dt * inner(div(u), q) *dx

b_interior = (
    inner(avg(grad(p_prev - p_trial)), jump(q, n)) * dS
    - inner(avg(grad(q)), jump(p_prev - p_trial, n)) * dS
    + beta/h * inner(jump(p_prev - p_trial, n), jump(q, n)) * dS
)

b_boundary = (inner(grad(p_prev - p_trial), q * n) * ds(2) + inner((p_prev - p_trial) * n, grad(q)) * ds(2) + beta_2/h * (p_prev - p_trial) * q * ds(2)
              + inner(grad(p_prev - p_trial), q * n) * ds(3)
    + inner((p_prev - p_trial) * n, grad(q)) * ds(3) + beta_2/h * (p_prev - p_trial) * q * ds(3) 
)

a2 = b_volume + b_boundary + b_interior
l2 =  0
function = a2 - l2

# Error with the below statement
solve(function == 0, p)
2

l2 cant be zero.

Secondly, you are setting up a non-linear problem in

but treating it as a linear problem as you use a TrialFunction in your variational form.

For further questions/clarifications, please add the full error message to your post as well.

3

Please find the code below along with the error message. I am trying to solve a linear problem using u_trial, but I am unable to understand why the error message ‘UFLException: Discontinuous type Jacobian must be restricted’ is appearing.
Code:

python
from fenics import *
import numpy as np

channel_length, channel_height = 1, 0.5
mesh = RectangleMesh(Point(0, 0), Point(channel_length, channel_height), 100, 100)

h = Constant(mesh.hmax())
n = FacetNormal(mesh)

V = VectorFunctionSpace(mesh, "DG", 2)
Q = FunctionSpace(mesh, "DG", 1)

u_trial = TrialFunction(V)
v = TestFunction(V)
u_prev = Function(V)
u_next  = Function(V)
u = Function(V)
p_prev = Function(Q)

Re = 0.01
e_field_ext = Constant((1e-3, 0))
k_electric = 1e-4
gamma = 500
gamma_2 = 1000

class InletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - 0) < DOLFIN_EPS and on_boundary

class OutletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - channel_length) < DOLFIN_EPS and on_boundary

inlet = InletBoundary()
outlet = OutletBoundary()

class DirichletBoundary_other(SubDomain):
    def inside(self, x, on_boundary):
        return (
            on_boundary
            and not (
                inlet.inside(x, on_boundary) or
                outlet.inside(x, on_boundary)
            )
        )

other = DirichletBoundary_other()

dx = Measure('dx', domain=mesh)
ds = Measure('ds', domain=mesh)
dS = Measure('dS', domain=mesh)

boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_markers.set_all(0)
other.mark(boundary_markers, 1)
inlet.mark(boundary_markers, 2)
outlet.mark(boundary_markers, 3)

ds = Measure("ds", domain=mesh, subdomain_data=boundary_markers)

dt = 0.001

#problem
a_volume = Re*(1/dt)*inner(u_trial - u_prev, v)*dx + dot(grad(p_prev),v)*dx + inner(grad(u_trial), grad(v))*dx 
a_interior = (
    - inner(avg(grad(u_trial)), outer(jump(v), n)) * dS
    - inner(avg(grad(v)), outer(jump(u_trial), n)) * dS
    + gamma/h * inner(jump(u_trial, n), jump(v, n)) * dS
)
a_boundary = (- inner(grad(u_trial), outer(v, n)) * ds(1) + inner(outer(u_trial, n), grad(v)) * ds(1) + gamma_2/h * dot(u_trial, v) * ds(1)
             )
function1 = a_volume + a_boundary + a_interior
a1 = lhs(function1)
l1 = rhs(function1)

#Error term
solve(a1 == l1, u)

Error:

python
 Calling FFC just-in-time (JIT) compiler, this may take some time.
Level 25:FFC:  Calling FFC just-in-time (JIT) compiler, this may take some time.
INFO:FFC:  Compiling form ffc_form_7692d3cd98a8ef442e8b0e28fa9c32b403118228
  
INFO:FFC:  Compiler stage 1: Analyzing form(s)
INFO:FFC:  -----------------------------------
DEBUG:FFC:    Preprocessing form using 'uflacs' representation family.
Discontinuous type Jacobian must be restricted.
ERROR:UFL_LEGACY:Discontinuous type Jacobian must be restricted.
---------------------------------------------------------------------------
UFLException                              Traceback (most recent call last)
<ipython-input-7-d0175a8e699d> in <cell line: 77>()
     75 
     76 #Error term
---> 77 solve(a1 == l1, u)

25 frames
/usr/local/lib/python3.10/dist-packages/ufl_legacy/log.py in error(self, *message)
    156         "Write error message and raise an exception."
    157         self._log.error(*message)
--> 158         raise self._exception_type(self._format_raw(*message))
    159 
    160     def begin(self, *message):

UFLException: Discontinuous type Jacobian must be restricted.
4

Unfortunately I can’t run your example on my current machine. But it looks like in your formulation for a there’s a restriction missing. It might be the variable h here,

    + gamma/h * inner(jump(u_trial, n), jump(v, n)) * dS

You could try something like h("+").

See also 1d dG advection/diffusion probem

5

Thank you very much for your reply. The code is now working with the following corrections:

a_interior = (                                                                 
    - inner(avg(grad(u_trial)), outer(jump(v), n("+"))) * dS
    - inner(avg(grad(v)), outer(jump(u_trial), n("+"))) * dS
    + gamma/h("+") * inner(jump(u_trial), jump(v)) * dS
)