Hey there community. I’m so new at FenicsX and I’m trying to solve the Taylor-Green Vortex problem. The velocity converges really well, but the pressure converges to something different and I don’t know how to correct it.
I’m using a monolithic schema of a IMEX Galerkin method. I left the code for everyone who wants or could help me to figure it out, or already has had the same problem as me 
.
In the picture attached I show on the left side the computed pressure field and in the right side the exact solution for the pressure field, both for the first time step (1e-5). I’m using the last fenics-dolfinx version (0.8.0) installed by conda-forge.
Thanks in advance.
#Taylor-Green Vortex
import numpy as np
import ufl
from dolfinx import fem, mesh, io
from mpi4py import MPI
from dolfinx import *
import dolfinx
from petsc4py import PETSc
from pathlib import Path
from dolfinx.mesh import create_unit_square
from basix.ufl import element, mixed_element
from dolfinx.fem.petsc import LinearProblem, assemble_matrix, assemble_vector, apply_lifting, create_vector, create_matrix,set_bc
from dolfinx.fem import Constant, Function, functionspace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical
from ufl import Measure, SpatialCoordinate, TestFunction, TrialFunction, div, exp, inner,dx,rhs,lhs,inner,dot,nabla_grad, TestFunctions, TrialFunctions
# Define the time-dependent function F_t
#define the exact solution for velocity
class vel_exac:
def __init__(self, t):
self.t = t
def eval(self, x):
values = np.vstack((np.sin(np.pi * x[0]) * np.cos(np.pi * x[1]) * F_t(t),
-np.sin(np.pi * x[1]) * np.cos(np.pi * x[0]) * F_t(t)),
dtype=PETSc.ScalarType)
return values
#define the exact solution for pressure
class press_exac:
def __init__(self, t):
self.t = t
def eval(self, x):
values = -0.5 * ((np.sin(2 * np.pi * x[0])**2) + (np.sin(2 * np.pi * x[1])**2)) * F_t(t)**2
return values
# Define the boundary condition
class Vel_BoundaryCondition:
def __init__(self, t):
self.t = t
def eval(self, x):
values = np.vstack((np.sin(np.pi * x[0]) * np.cos(np.pi * x[1]) * F_t(t),
-np.sin(np.pi * x[1]) * np.cos(np.pi * x[0]) * F_t(t)),
dtype=PETSc.ScalarType)
return values
def F_t(t):
return np.exp(-2 * np.pi**2 * t) # Example time-dependent function
def comp_and_save_exact_sol(t):
# Update the exact solution and saving
u_e.t = t; p_e.t = t
w_ex.sub(0).interpolate(u_e.eval) ; w_ex.sub(1).interpolate(p_e.eval)
u_ex_t = w_ex.sub(0).collapse() ; p_ex_t = w_ex.sub(1).collapse()
pressure_integral = assemble_scalar(form(p_ex_t * dx))
p_ex_t.vector[:] -= pressure_integral
u_ex.x.array[:] = u_ex_t.x.array[:] ; p_ex.x.array[:] = p_ex_t.x.array[:]
uex.write(t); pex.write(t)
return u_ex, p_ex
def all_walls(x):
return np.logical_or(np.logical_or(np.isclose(x[0], 0), np.isclose(x[0], 1)),
np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], 1)))
#Def point
def point(x):
point = [0.0, 0.0]
return np.logical_and(np.isclose(x[0], point[0]), np.isclose(x[1], point[1]))
# Number of elements in each direction of the mesh
n_elem = 80
msh = create_unit_square(MPI.COMM_WORLD, n_elem, n_elem)
#Approximation space polynomial degree
deg = 2
P2 = element("Lagrange", msh.basix_cell(), deg, shape=(msh.geometry.dim,))
P1 = element("Lagrange", msh.basix_cell(), deg-1)
#V, Q = functionspace(msh, P2), functionspace(msh, P1)
#defining mixed function space
W_el = mixed_element([P2, P1])
W = functionspace(msh, W_el)
# Create functions for the exact solutions
w_ex = Function(W) #; u_ex, p_ex = w_ex.split()
t = 0.0
u_e = vel_exac(t); p_e = press_exac(t)
w_ex.sub(0).interpolate(u_e.eval)
w_ex.sub(1).interpolate(p_e.eval)
u_ex = w_ex.sub(0).collapse() ; p_ex = w_ex.sub(1).collapse()
u_ex.name = "Velocity_exact"; p_ex.name = "Pressure_exact";
# Save the solution in BP format
output = Path("output")
output.mkdir(exist_ok=True)
#writing the exact solution in bp format
uex = io.VTXWriter(msh.comm, output / f"exact_solution_u.bp", u_ex)
pex = io.VTXWriter(msh.comm, output / f"exact_solution_p.bp", p_ex)
# integrate the pressure over the domain to subtract it to the pressure solution
pressure_integral = assemble_scalar(form(p_ex * dx))
pressure_integral_function = Constant(msh, PETSc.ScalarType(pressure_integral))
p_ex.vector[:] -= pressure_integral_function.value
uex.write(t); pex.write(t)
V,_ = W.sub(0).collapse()
Q,_ = W.sub(1).collapse()
#Facets and dofs walls
fdim = msh.topology.dim - 1
#Velocity BC on walls
facets_walls = dolfinx.mesh.locate_entities_boundary(msh, fdim, all_walls)
dofs_walls = dolfinx.fem.locate_dofs_topological((W.sub(0), V), 1, facets_walls)
u_bc_fes = Function(V)
u_values = vel_exac(t)
u_bc_fes.interpolate(u_values.eval)
bc_u = dirichletbc(u_bc_fes, dofs_walls, W.sub(0))
P_val = 0.0 #pressure Value at the point
point_located = dolfinx.mesh.locate_entities_boundary(msh, 0, point)
point_dofs = dolfinx.fem.locate_dofs_topological((W.sub(1),Q), 0, point_located)
press_fes = Function(Q)
press_fes.vector.set(P_val)
bc_p = dirichletbc(press_fes, point_dofs, W.sub(1))
bcs_up = [bc_u, bc_p]
#Initial condition
w_n = Function(W)
w_n.sub(0).interpolate(u_values.eval) #Exact solution at t=0
u_n, p_n = w_n.split()
# Save the solution in BP format
u_n_w = u_n.collapse(); p_n_w = p_n.collapse()
u_n_w.name = "Velocity"; p_n_w.name = "Pressure"
u_sol = io.VTXWriter(msh.comm, output / f"results_u.bp", u_n_w)
p_sol = io.VTXWriter(msh.comm, output / f"results_p.bp", p_n_w)
u_sol.write(t); p_sol.write(t)
# Save the boundary conditions in BP format
w_BC = Function(W)
w_BC.sub(1).interpolate(lambda x: x[0]*0 + 1) #Only to see the pressure BC
set_bc(w_BC.vector, bcs_up)
u_BC = w_BC.sub(0).collapse(); p_BC = w_BC.sub(1).collapse()
u_BC.name = "Velocity_BC"; p_BC.name = "Pressure_BC"
bc_u_w = io.VTXWriter(msh.comm,output / "BC_checking_u.bp", u_BC)
bc_p_w = io.VTXWriter(msh.comm,output / "BC_checking_p.bp", p_BC)
bc_u_w.write(t); bc_p_w.write(t)
#defininf Trial and Test functions
(u,p) = TrialFunctions(W)
(v,q) = TestFunctions(W)
dt = 1e-4
dx = Measure("dx",msh)
f = Constant(msh, PETSc.ScalarType((0, 0)))
nu = Constant(msh,PETSc.ScalarType(1))
tau = Constant(msh,PETSc.ScalarType(1/dt))
# Define the variational form
F1 = tau * dot(u,v)* dx # Mass term unknown
F1 += nu * inner(nabla_grad(u), nabla_grad(v)) * dx # Diffusive term
F1 += dot(dot(nabla_grad(u),u_n) , v) * dx # Convective term
F1 += dot(div(u),q) * dx - dot(p,div(v)) * dx # Continuity equation
F1 -= dot(f,v) * dx # RHS
F1 -= tau * dot(u_n,v)* dx # Mass term known
# Define the variational form for the LHS and RHS
a1 = form(lhs(F1)); l1 = form(rhs(F1))
#Assembling
#A = assemble_matrix(a1, bcs=bcs_up)
A = create_matrix(a1)
b = create_vector(l1)
# Create LU solver
solver = PETSc.KSP().create(msh.comm)
solver.setOperators(A)
solver.setType("preonly")
solver.getPC().setType("lu")
solver.getPC().setFactorSolverType("superlu_dist")
w_n1 = Function(W); u_n1, p_n1 = w_n1.split()
# Time-stepping
t = 0.0; t_end = 1
cnt = 0
# Define the error
L2_error_u = form(dot(u_n1 - u_ex, u_n1 - u_ex) * dx) #L2 error for the velocity
L2_error_p = form((p_n1 - p_ex)**2 * dx) #L2 error for the pressure
L2_history_w = [];
while t <= t_end:
t += dt; cnt += 1
t = round(t, 6)
print(f"t = {t:.5f}", end="\r")
# Update the exact solution and saving
u_ex, p_ex = comp_and_save_exact_sol(t)
# Update the time-dependent boundary condition
u_values.t = t; u_bc_fes.interpolate(u_values.eval)
bc_u = dirichletbc(u_bc_fes, dofs_walls, W.sub(0))
bc_p = dirichletbc(press_fes, point_dofs, W.sub(1))
bcs_up = [bc_u, bc_p]
# Assemble the matrix
A.zeroEntries()
assemble_matrix(A, a1, bcs=bcs_up)
A.assemble()
with b.localForm() as loc:
loc.set(0)
# Compute solution
assemble_vector(b, l1)
apply_lifting(b, [a1], [bcs_up])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, bcs_up)
solver.solve(b, w_n1.vector)
w_n1.x.scatter_forward()
# Split the solution and save it
u_n1, p_n1 = w_n1.split()
#compute the L2 norm of the error for u and p
error_L2_u = np.sqrt(msh.comm.allreduce(assemble_scalar(L2_error_u), op=MPI.SUM))
error_L2_p = np.sqrt(msh.comm.allreduce(assemble_scalar(L2_error_p), op=MPI.SUM))
print(f"Time: {t:.2e}; L2 error_u: {error_L2_u:.2e}; L2 error_p: {error_L2_p:.2e}")
L2_history_w.append([t, error_L2_u, error_L2_p])
#Update the solution
u_n.x.array[:] = u_n1.x.array[:]
p_n.x.array[:] = p_n1.x.array[:]
# Save the solution in BP format
u_temp = u_n.collapse(); p_temp = p_n.collapse()
u_n_w.x.array[:] = u_temp.x.array[:]; p_n_w.x.array[:] = p_temp.x.array[:]
u_sol.write(t); p_sol.write(t)
# Save the boundary conditions in BP format
set_bc(w_BC.vector, bcs_up)
u_BC.x.array[:] = w_BC.sub(0).collapse().x.array[:]
p_BC.x.array[:] = w_BC.sub(1).collapse().x.array[:]
bc_u_w.write(t); bc_p_w.write(t)
#only for solution comparation
if cnt == 10:
break
stop()
# Save the L2 error history
with open(output / "L2_error_history.txt", "w") as f:
f.write("Time L2_error_u L2_error_p\n")
for L2_error in L2_history_w:
f.write(f"{L2_error[0]} {L2_error[1]} {L2_error[2]}\n")
print('Computation completed ...')
uex.close(); pex.close()
u_sol.close(); p_sol.close()
bc_u_w.close(); bc_p_w.close()
solver.destroy()
b.destroy()
print('Solution written in output folder')