Hi everyone!
I want to use an RT0 variable to approximate the gradient of a DGP0 variable obtained from the mass-conservation equation.
The governing equations of \phi(\boldsymbol{x},t) and \boldsymbol{r}(\boldsymbol{x},t) are
\partial_{t} \phi+\nabla\cdot(\phi \boldsymbol{u})=0 in \Omega
\boldsymbol{r}=\nabla \phi in \Omega.
with boundary condition \nabla \phi\cdot \boldsymbol{n}=0 on \partial\Omega. Here \boldsymbol{u}(\boldsymbol{x},t) is the known velocity function used in the upwind operator.
We choose \phi\in DGP_{0} and \boldsymbol{r}\in RT_{0}, and use the upwind strategy: for any \psi\in DGP_{0} and \boldsymbol{s}\in RT_{0}, we have
step1: (\phi^{k+1}/\delta t, \psi)_{K}+(\phi^{k+1}_{upwind}(\boldsymbol{u}^{k}\cdot\boldsymbol{n}_{\partial K}), \psi)_{\partial K}=(\phi^{k}/\delta t, \psi)_{K}+(f_{source}^{k+1}, \psi)_{K},
where \phi^{k+1}_{upwind}=\phi^{k+1}|_{K^{+}} when \boldsymbol{u}^{k}\cdot\boldsymbol{n}_{\partial K}\geq 0, \phi^{k+1}_{upwind}=\phi^{k+1}|_{K^{-}} when \boldsymbol{u}^{k}\cdot\boldsymbol{n}_{\partial K}<0
step2: (\boldsymbol{r}^{k+1},\boldsymbol{s})_{K}=-(\phi^{k+1},\nabla\cdot\boldsymbol{s})_{K}+(\{\phi^{k+1} \boldsymbol{n}_{\partial K}\},[\boldsymbol{s}])_{\partial K}. with \boldsymbol{r}^{k+1}\cdot\boldsymbol{n}=0 on \partial\Omega.
The convergence order of \phi seems correct, but that of \boldsymbol{r} seems incorrect.
phi errors:
h | phi_L2 order | phi_Linf order | phi_H1 order
--------------------------------------------------------------------------------
1/4 | 9.5643e-02 | 2.6140e-01 | 9.9982e-01
1/8 | 3.9099e-02 1.2905 | 1.3968e-01 0.9041 | 9.9966e-01 0.0002
1/16 | 1.7238e-02 1.1815 | 6.8275e-02 1.0327 | 9.9965e-01 0.0000
1/32 | 8.0085e-03 1.1060 | 3.3997e-02 1.0059 | 9.9965e-01 0.0000
r errors:
h | r_L2 order | r_Linf order | div r L2 order
--------------------------------------------------------------------------------
1/4 | 5.6128e-01 | 1.2313e+00 | 6.5194e+00
1/8 | 4.6141e-01 0.2827 | 1.6404e+00 -0.4138 | 1.6878e+01 -1.3723
1/16 | 4.1555e-01 0.1510 | 2.0346e+00 -0.3108 | 3.6188e+01 -1.1004
1/32 | 3.9387e-01 0.0773 | 1.9944e+00 0.0288 | 7.3077e+01 -1.0139
At first, I suspected that the problem was caused by an incorrect implementation of \boldsymbol{r}.
Interestingly, when I replace the update of \phi with the interpolated exact solution, the computation of \boldsymbol{r} becomes correct.
"""update phi by exact solution"""
# phi_now.interpolate(lambda x: phi_e_expr(x, t_now))
r errors:
h | r_L2 order | r_Linf order | div r L2 order
--------------------------------------------------------------------------------
1/4 | 2.2790e-01 | 5.6240e-01 | 1.3033e+00
1/8 | 1.1347e-01 1.0061 | 2.7839e-01 1.0145 | 6.2398e-01 1.0625
1/16 | 5.6680e-02 1.0014 | 1.3889e-01 1.0031 | 3.0516e-01 1.0320
1/32 | 2.8333e-02 1.0003 | 6.9408e-02 1.0008 | 1.5093e-01 1.0157
This might be a finite element theory issue, or it could be a problem with my implementation.
I would appreciate any insights or explanations!
Below is my code.
import numpy as np
import sympy as sp
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import mesh
from dolfinx.fem import (Function, functionspace, dirichletbc, assemble_scalar,
form, locate_dofs_topological, Constant)
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, create_vector, apply_lifting, set_bc
from ufl import (CellDiameter, FacetNormal, TestFunction, TrialFunction, avg, jump,
nabla_grad, nabla_div, conditional, dot, dS, ds, dx, gt, inner, outer)
if np.issubdtype(PETSc.ScalarType, np.complexfloating):
print("Demo should only be executed with DOLFINx real mode")
exit(0)
def norm_l2(comm, v):
return np.sqrt(comm.allreduce(assemble_scalar(form(inner(v, v) * dx)), op=MPI.SUM))
def norm_h1(comm, v):
return np.sqrt(comm.allreduce(assemble_scalar(form(inner(nabla_grad(v), nabla_grad(v)) * dx)), op=MPI.SUM))
def norm_linf_phi(comm, u_h, u_exact_func, t_final, high_degree=5):
msh = u_h.function_space.mesh
V_high = functionspace(msh, ("Lagrange", high_degree))
u_h_high = Function(V_high)
u_h_high.interpolate(u_h)
u_exact_high = Function(V_high)
u_exact_high.interpolate(lambda x: u_exact_func(x, t_final))
local_max = np.max(np.abs(u_h_high.x.array - u_exact_high.x.array))
return comm.allreduce(local_max, op=MPI.MAX)
def norm_linf_r(comm, u_h, u_exact_func, t_final, high_degree=5):
msh = u_h.function_space.mesh
gdim = msh.geometry.dim
V_high = functionspace(msh, ("Lagrange", high_degree, (gdim,)))
u_h_high = Function(V_high)
u_h_high.interpolate(u_h)
u_exact_high = Function(V_high)
u_exact_high.interpolate(lambda x: u_exact_func(x, t_final))
local_max = np.max(np.abs(u_h_high.x.array - u_exact_high.x.array))
return comm.allreduce(local_max, op=MPI.MAX)
# TODO: exact solution(sympy)
sp_x, sp_y = sp.symbols('x y')
sp_t = sp.symbols('t')
"""variable"""
phi_sp = sp.cos(sp.pi * sp_t) * sp.cos(sp.pi * sp_x) * sp.cos(sp.pi * sp_y) * 0.45 + 0.5
u_sp = sp.Matrix([sp.cos(sp.pi * sp_t) * sp.sin(sp.pi * sp_x) * sp.sin(2 * sp.pi * sp_y),
-sp.cos(sp.pi * sp_t) * sp.sin(sp.pi * sp_y) * sp.sin(2 * sp.pi * sp_x)])
r_sp = sp.Matrix([sp.diff(phi_sp, sp_x),
sp.diff(phi_sp, sp_y)])
div_r_sp = sp.diff(r_sp[0], sp_x) + sp.diff(r_sp[1], sp_y)
_phi_e_expr = sp.lambdify((sp_x, sp_y, sp_t), phi_sp, 'numpy')
_u_e_expr = sp.lambdify((sp_x, sp_y, sp_t), u_sp, 'numpy')
_r_e_expr = sp.lambdify((sp_x, sp_y, sp_t), r_sp, 'numpy')
_div_r_e_expr = sp.lambdify((sp_x, sp_y, sp_t), div_r_sp, 'numpy')
"""mass eq"""
phi_u_sp = phi_sp * u_sp
rhs_phi_sp = sp.diff(phi_sp, sp_t) + sp.diff(phi_u_sp[0], sp_x) + sp.diff(phi_u_sp[1], sp_y)
_f_phi_expr = sp.lambdify((sp_x, sp_y, sp_t), rhs_phi_sp, 'numpy')
def phi_e_expr(x, t):
return _phi_e_expr(x[0], x[1], t)
def r_e_expr(x, t):
vals = _r_e_expr(x[0], x[1], t)
return np.array(vals).reshape(2, -1) if len(x.shape) > 1 else np.array(vals).flatten()
def div_r_e_expr(x, t):
return _div_r_e_expr(x[0], x[1], t)
def u_e_expr(x, t):
vals = _u_e_expr(x[0], x[1], t)
return np.array(vals).reshape(2, -1) if len(x.shape) > 1 else np.array(vals).flatten()
def f_phi_expr(x, t):
return _f_phi_expr(x[0], x[1], t)
def boundary_marker(x):
return np.isclose(x[0], 0.0) | np.isclose(x[0], 1.0) | np.isclose(x[1], 0.0) | np.isclose(x[1], 1.0)
left, right = 0, 1
bottom, top = 0, 1
t_final = 1
degree_phi = 0
degree_u = 1
mesh_i = np.array([4, 8, 16, 32])
mesh_j = mesh_i
mesh_t = (mesh_i**2).astype(int)
max_iter = np.size(mesh_i)
error_phi = np.zeros([max_iter, 3])
error_r = np.zeros([max_iter, 3])
for it in range(max_iter):
msh = mesh.create_rectangle(
comm=MPI.COMM_WORLD,
points=((left, bottom), (right, top)),
n=(mesh_i[it], mesh_j[it]),
cell_type=mesh.CellType.triangle,
)
# TODO: finite element space
P = functionspace(msh, ("Discontinuous Lagrange", degree_phi))
phi, psi = TrialFunction(P), TestFunction(P)
Q = functionspace(msh, ("Discontinuous Lagrange", degree_phi + 1))
K = functionspace(msh, ("Raviart-Thomas", degree_u))
r, s = TrialFunction(K), TestFunction(K)
# V = functionspace(msh, ("Brezzi-Douglas-Marini", degree_u))
V = functionspace(msh, ("Raviart-Thomas", degree_u + 1))
gdim = msh.geometry.dim
n = FacetNormal(msh)
nt = mesh_t[it]
dt = Constant(msh, t_final / nt)
t_now = 0.0
mass_err = np.zeros(nt)
momentum_err = np.zeros(nt)
fdim = msh.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(msh, fdim, boundary_marker)
r_D = Function(K)
r_D.interpolate(lambda x: r_e_expr(x, t_now))
boundary_vel_dofs_K = locate_dofs_topological(K, fdim, boundary_facets)
bc_r = dirichletbc(r_D, boundary_vel_dofs_K)
bcs_r = [bc_r]
r_now = Function(K)
f_phi = Function(Q)
f_phi.interpolate(lambda x: f_phi_expr(x, t_now))
phi_old = Function(P)
phi_now = Function(P)
phi_now.interpolate(lambda x: phi_e_expr(x, t_now))
u_old = Function(V)
u_now = Function(V)
u_now.interpolate(lambda x: u_e_expr(x, t_now))
"""lmbda = 1 (u_old dot n > 0); lmbda = 0 (u_old dot n <= 0)"""
lmbda = conditional(gt(dot(u_old, n), 0), 1, 0)
phi_uw = lmbda("+") * phi("+") + lmbda("-") * phi("-") # upwind
# TODO: mass equation
a_phi = form(
inner(phi / dt, psi) * dx
# convective terms
+ inner((dot(u_old, n))("+") * phi_uw, psi("+")) * dS
+ inner((dot(u_old, n))("-") * phi_uw, psi("-")) * dS
)
L_phi = form(
inner(phi_old / dt, psi) * dx
+ inner(f_phi, psi) * dx
)
# TODO: grad phi
a_r = form(
inner(r, s) * dx
)
L_r = form(
- inner(phi_now, nabla_div(s)) * dx
+ inner(avg(phi_now * n), jump(s)) * dS
)
A_phi = assemble_matrix(a_phi, [])
A_phi.assemble()
b_phi = create_vector(L_phi)
A_r = assemble_matrix(a_r, bcs_r)
A_r.assemble()
b_r = create_vector(L_r)
solver_phi = PETSc.KSP().create(msh.comm)
solver_phi.setOperators(A_phi)
solver_phi.setType(PETSc.KSP.Type.PREONLY)
solver_phi.getPC().setType(PETSc.PC.Type.LU)
solver_r = PETSc.KSP().create(msh.comm)
solver_r.setOperators(A_r)
solver_r.setType(PETSc.KSP.Type.PREONLY)
solver_r.getPC().setType(PETSc.PC.Type.LU)
phi_h = Function(P)
r_h = Function(K)
# TODO: time looping
for jt in range(nt):
t_now += dt.value
phi_old.x.array[:] = phi_now.x.array
phi_old.x.scatter_forward()
u_old.x.array[:] = u_now.x.array
u_old.x.scatter_forward()
# TODO: solve phi
f_phi.interpolate(lambda x: f_phi_expr(x, t_now))
A_phi.zeroEntries()
assemble_matrix(A_phi, a_phi, [])
A_phi.assemble()
with b_phi.localForm() as loc_b_phi:
loc_b_phi.set(0)
assemble_vector(b_phi, L_phi)
apply_lifting(b_phi, [a_phi], [[]])
b_phi.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
solver_phi.solve(b_phi, phi_h.x.petsc_vec)
phi_h.x.scatter_forward()
phi_now.x.array[:] = phi_h.x.array
phi_now.x.scatter_forward()
"""update phi by exact solution"""
# phi_now.interpolate(lambda x: phi_e_expr(x, t_now))
# TODO: solve grad phi
r_D.interpolate(lambda x: r_e_expr(x, t_now))
bc_r = dirichletbc(r_D, boundary_vel_dofs_K)
bcs_r = [bc_r]
A_r.zeroEntries()
assemble_matrix(A_r, a_r, bcs_r)
A_r.assemble()
with b_r.localForm() as loc_b_r:
loc_b_r.set(0)
assemble_vector(b_r, L_r)
apply_lifting(b_r, [a_r], [bcs_r])
b_r.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b_r, bcs_r)
solver_r.solve(b_r, r_h.x.petsc_vec)
r_h.x.scatter_forward()
r_now.x.array[:] = r_h.x.array
r_now.x.scatter_forward()
# TODO: update u
u_now.interpolate(lambda x: u_e_expr(x, t_now))
if msh.comm.rank == 0:
print('h = 1/%d, step = %d, t = %.4f' % (mesh_i[it], (jt + 1), t_now,))
P_e = functionspace(msh, ("Lagrange", degree_phi + 3))
V_e = functionspace(msh, ("Lagrange", degree_u + 3, (gdim,)))
phi_e = Function(P_e)
phi_e.interpolate(lambda x: phi_e_expr(x, t_now))
error_phi[it, 0] = norm_l2(msh.comm, phi_now - phi_e)
error_phi[it, 1] = norm_linf_phi(msh.comm, phi_now, phi_e_expr, t_now)
error_phi[it, 2] = norm_h1(msh.comm, phi_now - phi_e)
r_e = Function(V_e)
r_e.interpolate(lambda x: r_e_expr(x, t_now))
div_r_e = Function(P_e)
div_r_e.interpolate(lambda x: div_r_e_expr(x, t_now))
error_r[it, 0] = norm_l2(msh.comm, r_now - r_e)
error_r[it, 1] = norm_linf_r(msh.comm, r_now, r_e_expr, t_now)
error_r[it, 2] = norm_l2(msh.comm, nabla_div(r_now) - div_r_e)
if msh.comm.rank == 0:
print("phi errors:")
print(f"{'h':<8} | {'phi_L2':<13} {'order':<7} | {'phi_Linf':<13} {'order':<7} | {'phi_H1':<13} {'order':<7}")
print('-'*80)
for i in range(max_iter):
h_str = f"1/{mesh_i[i]}"
phi_l2, phi_linf, phi_h1 = error_phi[i, :]
if i == 0:
print(f"{h_str:<8} | {phi_l2:<13.4e} {'':<7} | {phi_linf:<13.4e} {'':<7} | {phi_h1:<13.4e} {'':<7}")
else:
phi_order_l2 = np.log(error_phi[i - 1, 0] / phi_l2) / np.log(2)
phi_order_linf = np.log(error_phi[i - 1, 1] / phi_linf) / np.log(2)
phi_order_h1 = np.log(error_phi[i - 1, 2] / phi_h1) / np.log(2)
print(f"{h_str:<8} | {phi_l2:<13.4e} {phi_order_l2:<7.4f} | {phi_linf:<13.4e} {phi_order_linf:<7.4f} | {phi_h1:<13.4e} {phi_order_h1:<7.4f}")
print("r errors:")
print(f"{'h':<8} | {'r_L2':<13} {'order':<7} | {'r_Linf':<13} {'order':<7} | {'div r L2':<13} {'order':<7}")
print('-'*80)
for i in range(max_iter):
h_str = f"1/{mesh_i[i]}"
r_l2, r_linf, div_r = error_r[i, :]
if i == 0:
print(f"{h_str:<8} | {r_l2:<13.4e} {'':<7} | {r_linf:<13.4e} {'':<7} | {div_r:<13.4e} {'':<7}")
else:
r_order_l2 = np.log(error_r[i - 1, 0] / r_l2) / np.log(2)
r_order_linf = np.log(error_r[i - 1, 1] / r_linf) / np.log(2)
div_r_order = np.log(error_r[i - 1, 2] / div_r) / np.log(2)
print(f"{h_str:<8} | {r_l2:<13.4e} {r_order_l2:<7.4f} | {r_linf:<13.4e} {r_order_linf:<7.4f} | {div_r:<13.4e} {div_r_order:<7.4f}")