I am attempting to set up a 1D (dG) 3D (FEM) advection-diffusion system and then run some convergence tests. I am running into what I believe may be some syntax errors. I have run multiple tests without my convection terms and the code works as a simple time-dependent diffusion problem. However when I add the advection term as shown in the code below, I get the error shown below the code chunk. Is there a particular syntax I must use when upwinding that I am not familiar with?
from math import gamma
from argparse import RawDescriptionHelpFormatter
from petsc4py import PETSc
from ast import Constant
from dolfin import *
from xii import *
import numpy as np
import time
kappa = 1.0
gamma_in = 1.0
gamma_out = 1.0
mesh_vec = [0,1,2,3]
T = 1.0
for i in mesh_vec:
ncell = 4*(2**i)
mesh_3d = BoxMesh(Point(-0.5,-0.5,-0.5),Point(0.5,0.5,0.5),ncell,ncell,ncell)
radius = 0.05
dt = 0.1/ncell
inv_dt = Constant(1.0/dt)
tpts = int(T/dt)
# === Perimeter and Area ===
D_perimeter = 2*np.pi*radius
D_area = np.pi*(radius**2)
# === Create Meshes ===
Omega = MeshFunction("size_t",mesh_3d,1,0)
CompiledSubDomain("near(x[0],0.0) && near(x[1],0.0)").mark(Omega,1)
Lambda = EmbeddedMesh(Omega,1)
# === Function Spaces ===
V3 = FunctionSpace(mesh_3d, "CG", 1)
V1 = FunctionSpace(Lambda, "DG", 1)
W = [V3, V1]
u3, u1 = map(TrialFunction, W)
v3, v1 = map(TestFunction, W)
# === Velocity terms ===
q1 = Constant(1.0)
q3 = as_vector((0,0,q1))
# === Initial concentration ===
cinlet = Constant(1.0)
# === Mesh Measures ===
dxOmega = Measure("dx", domain=mesh_3d)
dxLambda = Measure("dx", domain=Lambda)
dsLambda = Measure("ds", domain=Lambda)
n = FacetNormal(Lambda)
h = CellDiameter(Lambda)
alpha = Constant(30.0)
epsilon = Constant(-1.0)
perf1 = 1.0
perf3 = 1.0
# === Upwind ===
u1_up = u1('+')
v1_up = v1('+')
# === Cylinder ===
cylinder = Circle(radius=radius, degree=20)
# === Averages ===
u3_avg = Average(u3, Lambda, cylinder)
v3_avg = Average(v3, Lambda, cylinder)
# === Initialize time step ===
uh3d_prev = interpolate(Expression('0.0',degree=1), V3)
uh1d_prev = interpolate(Expression('0.0',degree=1), V1)
r = Expression('sqrt(pow(x[0] - 0.0, 2)+ pow(x[1] - 0.0,2))', degree = 7)
# === Begin time-stepping ===
t = 0.0
for step in range(tpts):
t += dt
print(f"[INFO] Time step {step+1}/{tpts}, t = {t:.3}")
# === Forcing terms ===
f3 = Expression('(r > radius ? area*(kappa/(1+kappa))*(1-radius*std::log(r/radius))*(sin(pi*x[2])+2) + (kappa/(1+kappa))*(1-radius*std::log(r/radius))*pow(pi,2)*sin(pi*x[2])*t :(kappa/(1+kappa))*pow(pi,2)*sin(pi*x[2]))*t + area*(kappa/(1+kappa))*(sin(pi*x[2])+2)',
degree = 4,
kappa = kappa,
radius = radius,
area = D_area,
r = r,
t = t)
f1 = Expression('(area*(sin(pi*x[2])+2) + area*sin(pi*x[2])*pow(pi,2)*t + (kappa/(1+kappa))*2*pi*radius*(sin(pi*x[2])+2))*t',
degree = 4,
radius = radius,
kappa = kappa,
area = D_area,
t = t)
# === Variational Forms ===
a00 = D_area * inv_dt * (inner(u3, v3) * dxOmega) + inner(grad(u3), grad(v3)) * dxOmega \
+ D_perimeter * kappa * inner(u3_avg, v3_avg) * dxLambda + inner(dot(q3,grad(u3)),v3) *dxOmega
a01 = -kappa * D_perimeter * inner(u1, v3_avg) * dxLambda
a10 = -kappa * D_perimeter * inner(u3_avg, v1) * dxLambda
a11 = D_area * inv_dt * inner(u1, v1) * dxLambda + D_area * inner(grad(u1), grad(v1)) * dxLambda \
- D_area * inner(avg(grad(u1)), jump(v1, n)) * dS \
- D_area * inner(avg(grad(v1)), jump(u1, n)) * dS \
+ (alpha/avg(h)) * inner(jump(u1), jump(v1)) * dS \
+ kappa * D_perimeter * inner(u1, v1) * dxLambda \
- D_area * q1 * inner(u1, grad(v1)) * dxLambda \
+ D_area * q1 * inner(u1_up, jump(v1, n)) * dS \
+ D_area * q1 * inner(v1_up, jump(u1, n)) * dS \
+ D_area * q1 * inner(u1,v1) * dsLambda(2) - D_area * inner(grad(u1),v1) * dsLambda(2) # boundary term
L0 = inner(f3, v3) * dxOmega + D_area * inv_dt * (inner(uh3d_prev, v3) * dxOmega)
L1 = inner(f1, v1) * dxLambda + D_area * inv_dt * (inner(uh1d_prev, v1) * dxLambda) \
+ D_area * q1 * inner(cinlet, v1) * dsLambda(1)
a = [[a00, a01], [a10, a11]]
L = [L0, L1]
bc3d = DirichletBC(V3, u3_exact, "on_boundary")
W_bcs = [bc3d, []] # No Dirichlet BCs on 1D
A, b = map(ii_assemble, (a, L))
A, b = apply_bc(A, b, W_bcs)
A, b = map(ii_convert, (A, b))
solver = PETScKrylovSolver()
solver.set_operators(A, A)
ksp = solver.ksp()
opts = PETSc.Options()
opts.setValue('ksp_type', 'gmres')
opts.setValue('ksp_norm_type', 'unpreconditioned')
opts.setValue('ksp_atol', 1E-14)
opts.setValue('ksp_rtol', 1E-30)
opts.setValue('ksp_monitor_true_residual', None)
opts.setValue('pc_type', 'hypre')
ksp.setFromOptions()
x = b.copy()
solver.solve(A,x,b)
wh = ii_Function(W)
wh.vector()[:] = x
uh3d, uh1d = wh
uh3d_prev.assign(uh3d)
uh1d_prev.assign(uh1d)
[INFO] Time step 1/40, t = 0.025
Shapes do not match: and .
ERROR:UFL_LEGACY:Shapes do not match: and .
UFLException Traceback (most recent call last)
/tmp/ipython-input-6-2479686055.py in <cell line: 0>()
154 + (alpha/avg(h)) * inner(jump(u1), jump(v1)) * dS
155 + kappa * D_perimeter * inner(u1, v1) * dxLambda
→ 156 - D_area * q1 * inner(u1, grad(v1)) * dxLambda
157 + D_area * q1 * inner(u1_up, jump(v1, n)) * dS
158 + D_area * q1 * inner(v1_up, jump(u1, n)) * dS \
2 frames
/usr/local/lib/python3.11/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: Shapes do not match: and .