Hello all,
I am currently in the process of trying to convert code for an incompressible flow mixed method from the Legacy version of Fenics to FenicsX. The legacy code is fully tested and works. However, I get completely different results when I run the FenicsX code. The legacy code with the correct results for time=2.0 s for the test case I am running is shown below:
from dolfin import *
from math import log
import numpy
parameters["std_out_all_processes"] = False;
parameters['ghost_mode'] = 'shared_facet'
########################
tol=1E-10
#######################
K=1 # Polynomial Order
################
dt= 0.1# Time Step
Time= 40 # Total Time
n_u = 0.001568 #0.00001568
rho = 1
########################
count = 0 # File output counter
###################### Read Mesh in and define boundaries
mesh = Mesh()
hdf = HDF5File(mesh.mpi_comm(), "Channel_mesh_fine.h5", "r")
hdf.read(mesh, "/mesh", False)
cd = MeshFunction("size_t", mesh,2)
hdf.read(cd, "/cd")
fd = MeshFunction("size_t", mesh,1)
hdf.read(fd, "/fd")
ds = Measure('ds', domain=mesh, subdomain_data=fd)
###################### Define Function Spaces and Elements
U=VectorFunctionSpace(mesh,'CG',K+1) # for velocity
P=FunctionSpace(mesh,'CG',K) # for pressure
U_el=VectorElement("CG", mesh.ufl_cell(), K+1) ### CG element is choosen for velocity
P_el= FiniteElement("CG", mesh.ufl_cell(), K) ### CG space is choosen for pressure
W_1=MixedElement([U_el, P_el])
W=FunctionSpace(mesh, W_1)
du=TrialFunction(W)
up=Function(W)
u,p=split(up)
vq=TestFunction(W)
v,q=split(vq)
####################### Initial Conditions
Velocity=Expression(("0.05", "0"),degree=1)
u_1=Function(U)
u_2=Function(U)
u_3=Function(U)
u_1=interpolate(Velocity,U)
u_2=interpolate(Velocity,U)
u_3=interpolate(Velocity,U)
############################################# BCs
# Velocity BCs
UNoslip=DirichletBC(W.sub(0),Constant((0.0, 0.0)),fd,3) # Wall no-slip
UNoslip2=DirichletBC(W.sub(0),Constant((0.0, 0.0)),fd,4) # Wall no-slip
UbcIn=DirichletBC(W.sub(0),Constant((0.05, 0.0)),fd,1)
# Pressure BCs
PbcOut = DirichletBC(W.sub(1), Constant(0.0), fd, 2)# kg/m*s^2
bcs=[UNoslip,UNoslip2, UbcIn, PbcOut]
#############################################
cons=6.0/11.0*dt # BDF3 hf Coefficient
iden=Identity(2) # Identity Matrix
########################
L_1=dot(u,v)*dx-18.0/11.0*dot(u_3,v)*dx+9.0/11.0*dot(u_2,v)*dx-2.0/11.0*dot(u_1,v)*dx
########################
L_2=-cons*div(u)*q*dx
#########################
L_3=cons*dot(dot(u,nabla_grad(u)),v)*dx\
+cons*1/2*dot(div(u)*u,v)*dx\
-cons*(1/rho)*p*div(v)*dx
###########################
L_4=cons*n_u*inner(nabla_grad(u)+nabla_grad(u).T-2/3*div(u)*iden,nabla_grad(v))*dx
##############################
L=L_1+L_2+L_3+L_4
t_i=2*dt
###################### Set Up Solver
J= derivative(L,up,du)
problem = NonlinearVariationalProblem(L,up,bcs,J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm["newton_solver"]["absolute_tolerance"] = tol
prm["newton_solver"]["relative_tolerance"] = tol
prm["newton_solver"]["maximum_iterations"] = 25
#prm["newton_solver"]["relaxation_parameter"] = 1.0
#prm["newton_solver"]["linear_solver"] = "bicgstab"
#prm["newton_solver"]["linear_solver"] = "petsc"
#prm["newton_solver"]["preconditioner"] = "icc"
prm["newton_solver"]["linear_solver"] = "mumps" #most stable
solver.parameters["newton_solver"]["error_on_nonconvergence"] = False
while(t_i<Time):
t_i+=dt
print(t_i)
solver.solve()
_u,_p=up.split(True)
u_1.assign(u_2)
u_2.assign(u_3)
u_3.assign(_u)
if(count % 3 == 0):
File("../scratch/Velocity_"+"Time_"+str(t_i)+"_.pvd") << _u
File("../scratch/Pressure_" +"Time_"+str(t_i)+"_.pvd") << _p
count = count + 1
quit()
However, when I run the FenicsX code shown which has the same exact weak equation as the code above and essentially follows the code above with the FenicsX syntax, I get completely different results. It almost looks like the Reynolds number is way higher for the flow below. The boundary conditions seem to be implemented correctly as the output implements the no slip, inlet, and outlet (not shown below but there is 0 pressure in the outlet in pressure output) correctly, and the initial conditions are also fine. The code along with the output after 2s is shown below:
from petsc4py import PETSc
from dolfinx.fem.petsc import set_bc
from math import log
import numpy
#from ufl import (TestFunction, TrialFunction, div, dot, dS, ds, dx, grad, nabla_grad, inner, outer, split, Identity)
from ufl import *
from basix.ufl import element, mixed_element
from dolfinx import fem, la, default_scalar_type, io, nls, log
from dolfinx.fem import (
Constant,
Expression,
Function,
dirichletbc,
DirichletBC,
extract_function_spaces,
form,
functionspace,
locate_dofs_topological,
)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.io import XDMFFile, VTXWriter
from mpi4py import MPI
from pathlib import Path
from dolfinx.io import gmshio
from petsc4py.PETSc import ScalarType
################# Set Key Solver Parameters
tol = 1e-10 # Tolerance
K = 1 # Polynomial Degree
dt = 0.1 # Time Step
Time = 30 # Simulation time
count = 1; # Time loop counter
################# Set Key Solver Parameters
n_u = 0.001568 # Kinematic Viscosity (extra high to emphasize viscous effects)
rho = 1 # Approx Density of air
vel_i = 0.05; # Initial Velocity (set equal to velocity over the domain)
vel_in = 0.05 # Inlet Velocity
press_out = 0.0; # Outlet Pressure
################# Read in Mesh
with XDMFFile(MPI.COMM_WORLD, "Channel_mesh_fine.xdmf", "r") as xdmf:
mesh = xdmf.read_mesh(name="Grid")
cd = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)
with XDMFFile(MPI.COMM_WORLD, "Channel_mesh_fine_facet.xdmf", "r") as xdmf:
fd = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
################# Define Finite Elements
U_el = element("CG", mesh.basix_cell(), K + 1, shape = (mesh.geometry.dim,)) # Velocity Element
P_el = element("CG", mesh.basix_cell(), K) # Pressure Element
################# Define Function Spaces (For initial Conditions)
U = functionspace(mesh, U_el) # Velocity Space
P = functionspace(mesh, P_el) # Pressure Space
W_el = mixed_element([U_el, P_el]) # Define Mixed element
W = functionspace(mesh, W_el) # Define Mixed Function Space
################# Define Trial and Test Functions
du = TrialFunction(W) # For Jacobian
up = Function(W)
u,p = split(up)
vq = TestFunction(W)
v,q = split(vq)
################# Define Function Spaces for Runge-Kutta Terms
u_1 = Function(U)
u_2 = Function(U)
u_3 = Function(U)
################ Set up Initial Conditions using initial value for each RK Term
for i, vel in enumerate([vel_i, 0.0]):
u_1.sub(i).interpolate(lambda x, vel_in= vel: numpy.full(x.shape[1], vel_in, dtype=numpy.float64))
u_2.sub(i).interpolate(lambda x, vel_in = vel: numpy.full(x.shape[1], vel_in, dtype=numpy.float64))
u_3.sub(i).interpolate(lambda x, vel_in = vel: numpy.full(x.shape[1], vel_in, dtype=numpy.float64))
#### Check if Velocity was assigned properly
folder = Path("results")
folder.mkdir(exist_ok=True, parents=True)
vtx_u = VTXWriter(mesh.comm, "Vel_check.bp", [u_1], engine="BP4")
vtx_u.write(0.0)
################ Define Boundary Conditions
# Legacy BCs
#UNoslip=DirichletBC(W.sub(0),Constant((0.0, 0.0)),fd,3) # Wall no-slip
#UNoslip2=DirichletBC(W.sub(0),Constant((0.0, 0.0)),fd,4) # Wall no-slip
#UbcIn=DirichletBC(W.sub(0),Constant((0.05, 0.0)),fd,1) # Inlet Velocity
#PbcOut = DirichletBC(W.sub(1), Constant(0.0), fd, 2)# Outlet Pressure
#BCs=[UNoslip,UNoslip2, UbcIn, PbcOut]
# DolfinX BCs
tdim = mesh.topology.dim
fdim = tdim - 1
def inflow(x):
return numpy.stack((numpy.full(x.shape[1], 0.05), numpy.zeros(x.shape[1])))
def noslip(x):
return numpy.stack((numpy.zeros(x.shape[1]), numpy.zeros(x.shape[1])))
def outlet(x):
return numpy.zeros(x.shape[1])
ID_Wall_1 = 3
ID_Wall_2 = 4
ID_In = 1
ID_Out = 2
V_sub, _ = W.sub(0).collapse()
wall_1_dofs = locate_dofs_topological((W.sub(0),V_sub), fdim, fd.find(ID_Wall_1))
wall_2_dofs = locate_dofs_topological((W.sub(0),V_sub), fdim, fd.find(ID_Wall_2))
in_dofs = locate_dofs_topological((W.sub(0),V_sub), fdim, fd.find(ID_In))
out_dofs = locate_dofs_topological((W.sub(1)), fdim, fd.find(ID_Out))
u_wall_1 = Function(V_sub)
u_wall_2 = Function(V_sub)
u_wall_1.interpolate(noslip)
u_wall_2.interpolate(noslip)
bc_wall_1 = dirichletbc(u_wall_1, wall_1_dofs, W.sub(0))
bc_wall_2 = dirichletbc(u_wall_2, wall_2_dofs, W.sub(0))
u_in = Function(V_sub)
u_in.interpolate(inflow)
bc_inlet = dirichletbc(u_in, in_dofs, W.sub(0))
bc_outlet = dirichletbc(ScalarType(press_out), out_dofs, W.sub(1))
BCs = [bc_wall_1, bc_wall_2, bc_inlet, bc_outlet]
################ Weak Formulation Constants
cons=6.0/11.0*dt # BDF3 derivative coefficient
iden = Identity(2) # Identity Matrix
################ Weak Formulation
L_1=dot(u,v)*dx-18.0/11.0*dot(u_3,v)*dx+9.0/11.0*dot(u_2,v)*dx-2.0/11.0*dot(u_1,v)*dx
L_2=-cons*div(u)*q*dx
L_3=cons*dot(dot(u,nabla_grad(u)),v)*dx\
+cons*1/2*dot(div(u)*u,v)*dx\
-cons*(1/rho)*p*div(v)*dx
L_4=cons*n_u*inner(nabla_grad(u)+nabla_grad(u).T-2/3*div(u)*iden,nabla_grad(v))*dx\
L=L_1+L_2+L_3+L_4
t_i=2*dt
################ Set Solver Definition
Jac = derivative(L,up,du)
problem = NonlinearProblem(L,up,bcs=BCs, J = Jac)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = tol
solver.atol = tol
solver.report = True
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.INFO)
################ Loop over Time
while(t_i < Time):
t_i += dt
print(t_i)
n, converged = solver.solve(up)
assert (converged)
print(f"Number of interations: {n:d}")
u_ = up.sub(0).collapse()
p_ = up.sub(1).collapse()
u_.x.scatter_forward()
p_.x.scatter_forward()
########## Update stored values
u_1.x.array[:] = u_2.x.array[:]
u_2.x.array[:] = u_3.x.array[:]
u_3.x.array[:] = u_.x.array[:]
########### Print out values
if(count % 5 == 0 or count < 4):
folder = Path("results")
folder.mkdir(exist_ok=True, parents=True)
vtx_u = VTXWriter(mesh.comm, "Velocity_"+"Time_"+str(t_i)+"_.bp", [u_], engine="BP4")
vtx_p = VTXWriter(mesh.comm, "Pressure_"+"Time_"+str(t_i)+"_.bp", [p_], engine="BP4")
vtx_p.write(t_i)
vtx_u.write(t_i)
count += 1
quit()
I am not sure what the error is, so I appreciate any and all help. It might be related to the syntax switch for the mixed element, the solver, or the solution update, but I am not sure.Thanks!