I am trying to solve the unsteady viscous 2D Burgers equation with periodic BCs in the x-direction and Dirichlet BCs in the y-direction in a space domain [0,pi]x[0,2] and a time domain [0,1] with a time step dt=1e-4. I use the Crank-Nicholson scheme for time integration. The solution is defined in a VectorFunctionSpace and the initial condition and forcing term are defined from a DNS dataset such that the obtained FEniCS solution represents the spanwise and wall-normal velocity components from a Turbulent Channel Flow (TCF) at a Re_tau=180.
I am using linear triangular elements and thus I am expecting to have a convergence rate of O(h^2) where h is the element’s length. I use h-refinement with a factor of 2 across all spatial resolutions ranging from 32x16 to 512x256. However, that is not the case as the presented convergence is approximately 1.3. Any idea on why this is happening? I have performed a test case similar but in 1D from the same dataset (to just obtain the wall-normal component) and the convergence rate was very close to the expected.
Below is the implemented FEM solver for the 2D case:
from dolfin import *
import numpy as np
import os
class DNSDataExpression(UserExpression):
def __init__(self, ww, vv, **kwargs):
self.ww = ww
self.vv = vv
super().__init__(**kwargs)
def eval(self, values, x):
z_index = int((x[0] / L_z) * (self.ww.shape[1] - 1))
y_index = int((x[1] / L_y) * (self.ww.shape[0] - 1))
values[0] = self.ww[y_index, z_index]
values[1] = self.vv[y_index, z_index]
def value_shape(self):
return (2,)
class PeriodicBoundary(SubDomain):
def inside(self, x, on_boundary):
return bool(x[0] < DOLFIN_EPS and x[0] > -DOLFIN_EPS and on_boundary)
def map(self, x, y):
y[0] = x[0] - L_z
y[1] = x[1]
class DirichletBoundary(SubDomain):
def inside(self, x, on_boundary):
return bool((x[1] < DOLFIN_EPS or x[1] > (L_y - DOLFIN_EPS)) and on_boundary)
def InterpolateDNSData(N_zr, N_yr, f1, f2, dim_2D = True):
if dim_2D:
List item
f1 = f1[::N_yr, ::N_zr]
f2 = f2[::N_yr, ::N_zr]
else:
f1 = f1[:, ::N_yr, ::N_zr]
f2 = f2[:, ::N_yr, ::N_zr]
return f1, f2
def SolverParameters(solver):
# Parameters for non-linear solvers:
prm = solver.parameters
prm["nonlinear_solver"] = "snes"
prm["snes_solver"]["line_search"] = "bt"
prm["snes_solver"]["linear_solver"] = "mumps"
prm["snes_solver"]["krylov_solver"]["nonzero_initial_guess"] = True
prm["snes_solver"]["report"] = False
set_log_active(False)
return solver
N_z = 32 # Number of elements in spanwise z-direction
N_y = N_z // 2 # Number of elements in wall-normal y-direction
L_z = 3.14159 # Spanwise length
L_y = 2 # Wall-normal length
z = np.linspace(0, L_z, N_z + 1) # Coordinates in z-direction
y = np.linspace(0, L_y, N_y + 1) # Coordinates in y-direction
degree = 1 # Degree of the finite element basis functions
Re = 180 # Reynolds number
# Time integration parameters
T = 1.0
dt = 1e-4
t = np.linspace(0, T, int(T / dt) + 1)
N_t = len(t)
theta = 0.5
print('\nNumber of elements in z-direction: ', N_z)
print('Number of elements in y-direction: ', N_y)
N_zz = 512
N_yy = N_zz // 2
# Load DNS data
print('\nLoading DNS data')
filename = 'data/DNS/burgers_2D_forcing.npz'
file_path = os.path.join(os.path.dirname(__file__), filename)
data = np.load(file_path)
f_zz = data['f_z']
f_yy = data['f_y']
filename = 'data/DNS/burgers_2D_IC.npz'
file_path = os.path.join(os.path.dirname(__file__), filename)
data = np.load(file_path)
ww = data['w']
vv = data['v']
# Interpolation in space to have same resolution as FEM solver
N_zr = int(N_zz / N_z)
N_yr = int(N_yy / N_y)
ww, vv = InterpolateDNSData(N_zr, N_yr, ww, vv)
f_zz, f_yy = InterpolateDNSData(N_zr, N_yr, f_zz, f_yy, dim_2D = False)
# Create mesh and define function space
rectangle_mesh = RectangleMesh(Point(0, 0), Point(L_z, L_y), N_z, N_y)
V = VectorFunctionSpace(rectangle_mesh, 'CG', degree, constrained_domain = PeriodicBoundary())
# Define Dirichlet BC
u_D = Constant((0, 0))
bc = DirichletBC(V, u_D, DirichletBoundary())
bcs = [bc]
# Define necessary functions
u = Function(V)
u_0 = Function(V)
f = Function(V)
f_0 = Function(V)
v = TestFunction(V)
# Initial condition and forcing terms from DNS data
u_e = DNSDataExpression(ww, vv, degree = degree)
u.interpolate(u_e)
u_0.assign(u)
f_e = DNSDataExpression(f_zz[0], f_yy[0], degree = degree)
f.interpolate(f_e)
f_0.assign(f)
# Define variational problem
F_theta = dot(dot(u, nabla_grad(u)), v) + inner(grad(u), grad(v)) / Re - dot(f, v)
F_1_theta = dot(dot(u_0, nabla_grad(u_0)), v) + inner(grad(u_0), grad(v)) / Re - dot(f_0, v)
F = (dot((u - u_0) / dt, v)) * dx + theta * F_theta * dx + (1 - theta) * F_1_theta * dx
dF_du = derivative(F, u)
problem = NonlinearVariationalProblem(F, u, bcs, dF_du)
solver = NonlinearVariationalSolver(problem)
solver = SolverParameters(solver)
u_history = [u.copy(deepcopy=True)]
print('\nFEM solver started')
for i in range(1, N_t):
# Update source term
f_e = DNSDataExpression(f_zz[i], f_yy[i], degree = degree)
f.interpolate(f_e)
# Solve primal variational problem
solver.solve()
u_history.append(u.copy(deepcopy=True))
u_0.assign(u)
f_0.assign(f)