Hi, sure, I am importing a vtk mesh:
I tried the same code using a box generated by dolfinx and the code works. If I use the mesh I get this for the velocity magnitude:
The boundary conditions should be no slip at the cone and top cylinder (not the top cap) boundary facets, and at the bottom cap (which is not working at all)
Thank you for your time,
Alex
This is the running code:
vol_mesh = meshio.read('mesh.vtk')
points = vol_mesh.points
cells = vol_mesh.cells[1].data
c_el = ufl.Mesh(basix.ufl.element("Lagrange", "tetrahedron", 1, shape=(points.shape[1],)))
msh = mesh.create_mesh(MPI.COMM_WORLD, cells, points, c_el)
fdim = msh.topology.dim - 1
msh.topology.create_connectivity(fdim, fdim+1)
# Geometry variables
h_array_1 = 0.02
h_array_2 = 0.035
h_gap = 0.03
rho_air = fem.Constant(msh, 1.225)
mu_air = fem.Constant(msh, 1.8e-3)
# Function to mark noslip
def noslip_boundary1(x):
return np.logical_and(x[2] < h_array_2, x[2] > -h_array_1)
# Function to mark noslip
def noslip_boundary2(x):
return np.isclose(x[2], (h_array_2 + h_gap))
# the same as old ufl.VectorElement
# the same as old ufl.FiniteElement
P2 = basix.ufl.element("Lagrange", msh.topology.cell_name(), 2, shape=(msh.geometry.dim, ))
P1 = basix.ufl.element("Lagrange", msh.topology.cell_name(), 1)
V, Q = fem.functionspace(msh, P2), fem.functionspace(msh, P1)
# Create the function space
TH = basix.ufl.mixed_element([P2, P1])
W = fem.functionspace(msh, TH)
W0, _ = W.sub(0).collapse()
# No slip boundary condition 1
noslip = fem.Function(V)
facets1 = mesh.locate_entities_boundary(msh, fdim, noslip_boundary1)
dofs1 = fem.locate_dofs_topological((W.sub(0), V), fdim, facets1)
bc0 = fem.dirichletbc(noslip, dofs1, W.sub(0))
# No slip boundary condition 2
facets2 = mesh.locate_entities_boundary(msh, fdim, noslip_boundary2)
dofs2 = fem.locate_dofs_topological((W.sub(0), V), fdim, facets2)
bc1 = fem.dirichletbc(noslip, dofs2, W.sub(0))
# Since for this problem the pressure is only determined up to a constant, we pin it at a point
zero = fem.Function(Q)
zero.x.array[:] = 0
dofs = fem.locate_dofs_geometrical(
(W.sub(1), Q), lambda x: np.isclose(x.T, [0.0, 0.0, -h_array_1], atol=0.01).all(axis=1))
bc2 = fem.dirichletbc(zero, dofs, W.sub(1))
def body_force_fun(x):
F_max = 2.0e-2
z_offset = 0.02
r_xy = np.sqrt(x[0]**2 + x[1]**2)
a = 0.005
max_val = a * np.exp(-0.5)
Fr_val = np.exp(-r_xy**2/(2*a**2)) / max_val
b = 0.005
Fz_val = np.exp(-(x[2]-z_offset)**2/(2*b**2))
F_mag = F_max * Fr_val * Fz_val
# -ve sign to get clockwise rotation from array top
return -F_mag * np.vstack((x[1], -x[0], np.zeros_like(x[2])))
Force_per_unit_vol = fem.Function(V)
Force_per_unit_vol.interpolate(lambda x: body_force_fun(x))
# Collect Dirichlet boundary conditions
bcs = [bc0, bc1, bc2]
# Define variational problem
w = fem.Function(W)
(u, p) = ufl.split(w)
(v, q) = ufl.TestFunctions(W)
# Tentative velocity step
# https://github.com/bragostin/Fenics/blob/main/HTC_SquareDuct_Steady_Neumann_Minimal_Pub.py
# https://math.stackexchange.com/questions/1952588/
# weak-form-of-steady-navier-stokes-equations-with-special-boundary-condition
F = rho_air*inner(grad(u)*u, v)*dx + \
mu_air*inner(grad(u), grad(v))*dx + \
inner(Force_per_unit_vol, v)*dx - \
inner(p, div(v))*dx - \
div(u)*q*dx
dw = ufl.TrialFunction(W)
dF = ufl.derivative(F, w, dw)
problem = NonlinearProblem(F, w, bcs=bcs, J=dF)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
#solver.max_it = 2
solver.report = True
solver.error_on_nonconvergence = False
# Custom solver options (not needed for code readibility but makes code 9x faster)
opts = PETSc.Options() # type: ignore
ksp = solver.krylov_solver
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()
# Compute the solution
log.set_log_level(log.LogLevel.INFO)
start = time.time()
solver.solve(w)
end = time.time()
print('\n')
print(msh.topology.index_map(msh.topology.dim).size_local)
print((end - start) / 60)
# Split the mixed solution and collapse
u = w.sub(0).collapse()
p = w.sub(1).collapse()
# Write the solution to file
with io.XDMFFile(MPI.COMM_WORLD, "pressure.xdmf", "w") as pfile_xdmf:
p.x.scatter_forward()
pfile_xdmf.write_mesh(msh)
pfile_xdmf.write_function(p)
with io.XDMFFile(MPI.COMM_WORLD, "velocity.xdmf", "w") as ufile_xdmf:
u.x.scatter_forward()
P1 = basix.ufl.element("Lagrange", msh.basix_cell(), 1, shape=(msh.geometry.dim,))
u1 = fem.Function(fem.functionspace(msh, P1))
u1.interpolate(u)
ufile_xdmf.write_mesh(msh)
ufile_xdmf.write_function(u1)
with io.XDMFFile(MPI.COMM_WORLD, "body_force.xdmf", "w") as ufile_xdmf:
Force_per_unit_vol.x.scatter_forward()
P1 = basix.ufl.element("Lagrange", msh.basix_cell(), 1, shape=(msh.geometry.dim,))
u1 = fem.Function(fem.functionspace(msh, P1))
u1.interpolate(Force_per_unit_vol)
ufile_xdmf.write_mesh(msh)
ufile_xdmf.write_function(u1)
