1d dG advection/diffusion probem

dolfinx
1

Hi, I’m new to Fenics and dG implementation but trying to construct a 1D dG advection/diffusion code using a Robin bc at the inlet and Neumann at the outlet. We are assuming that there is a given velocity prescribed at the inlet and propagated to the outlet (from right to left). I keep getting the error “UFLException: Discontinuous type Argument must be restricted.” and do not understand what this means or what line it is pertaining to. I know it occurs when I try to use solve() but do not know where the issue lies within this command. I have included the code and error output below. Thanks in advance.

import copy
import datetime
import importlib
import json
import matplotlib.pyplot as plt
import meshio
import networkx as nx
import numpy as np
import scipy
import vtk

from dolfin import *
from graphnics import *
from matplotlib.gridspec import GridSpec
from mpl_toolkits.mplot3d import Axes3D
from scipy.spatial import cKDTree
from vtk.util.numpy_support import vtk_to_numpy
from xii import *

# Mesh

mesh = UnitIntervalMesh(50)

# Function spaces

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

# Veclocity: set as constant for now
u_constant = Constant(1.0)
u = u_constant

cn = Constant(1.0)
cinlet = Constant(1.0)

# Trial and test functions
c1 = TrialFunction(V_dg)
w  = TestFunction(V_dg)

# Diffusion coefficient
gamma = Constant(1.0)

# Time step
dt = Constant(0.01)

# Source term
c = cn/dt
g = Constant(1.0) 

# Penalty term
alpha = Constant(1.0)

# Mesh-related functions
n = FacetNormal(mesh)
h = CellDiameter(mesh)

# Upwind
uv = as_vector((u,))
un = (dot(uv,n) + abs(dot(uv,n)) )/2.0

# Define boundary subdomains
class Left(SubDomain):
  def inside(self, x, on_boundary):
            return on_boundary and near(x[0], 0)

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

left = Left()
right = Right()


# Mark boundaries
boundaries = MeshFunction('size_t', mesh,mesh.topology().dim()-1)
boundaries.set_all(0)
left.mark(boundaries,0)
right.mark(boundaries,1)

dxLambda = Measure('dx', domain = mesh)
dsLambda = Measure('ds', domain = mesh, subdomain_data = boundaries)

# Bilinear form

a_vol = dot(w,c1/dt)*dxLambda + dot(grad(w),grad(c1)-uv*c1)*dxLambda


a_fac = - dot(avg(grad(w)), jump(c1, n))*dS \
            - dot(jump(w, n), avg(grad(c1)))*dS \
            + (alpha('+')/h('+'))*dot(jump(w, n), jump(c1, n))*dS \
            + gamma('+')*dot(grad(w),grad(c1('+')))*dS


a_vel = dot(jump(w), un('+')*c1('+') - un('-')*c1('-'))*dS \
            +dot(w,un('-')*c1('-'))*dsLambda(1)

a     = a_vol + a_fac + a_vel

# Linear form

L = dot(w,c)*dxLambda + u*cinlet*w*dsLambda(0) + dot(g,w)*dxLambda

# Solution function

ch = Function(V_dg)


solve(a==L, ch)

Error message:

Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 98)
        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_8d578a6120e21c5936ce97e31284171d8c8980e6
        
INFO:FFC:        Compiler stage 1: Analyzing form(s)
INFO:FFC:        -----------------------------------
DEBUG:FFC:          Preprocessing form using 'uflacs' representation family.
Discontinuous type Argument must be restricted.
ERROR:UFL_LEGACY:Discontinuous type Argument must be restricted.
---------------------------------------------------------------------------
UFLException                              Traceback (most recent call last)
<ipython-input-7-706dc4a8007f> in <cell line: 125>()
    123 #solve(a == L, ch, bc, solver_parameters={"linear_solver": "mumps"})
    124 #solver = LUSolver(a==L, "mumps")
--> 125 solve(a==L, ch)
    126 

27 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 Argument must be restricted.
2

Unfortunately I don’t have a working legacy dolfin environment on my machine. But at first glance I notice:

            +dot(w,un('-')*c1('-'))*dsLambda(1)

UFL requires that coefficients be restricted to a particular side of the facet on interior facets (as you have done with ("+"), ("-"), jump(), avg(), etc.); however, on exterior facets no such restriction is required. So perhaps you could audit your formulation and make sure you adhere to this guideline. For example in the case above I would try:

            +dot(w,un*c1)*dsLambda(1)

If I can find a working docker image for legacy dolfin I’ll try your code and update this post.

3

Thanks to @dokken I managed to get a working image of legacy dolfin. The other offending line in your code is

+ gamma('+')*dot(grad(w),grad(c1('+')))*dS

You need to choose how to restrict w in this case.

4

Thank you so much for your reply! I was able to get it working with these fixes and I think now I have a better understanding of when to restrict vs when not to restrict. Thank you again!