Hello, According to
I wanted to ask if there is a more efficient way to implement the IPSCR code in the IPCSR benchmark. Although what I have done does not give me any errors, and what’s more, the IPCSR scheme obtained good results with an open boundary (u and p homogeneous Neumann conditions), which were not obtained with IPCS. However, when applied to an extensive domain (high Reynolds number), there were no improvements. Does anyone have any suggestions? SUPG and LSIC have already been implemented.
p_.x.petsc_vec.axpy(1, phi.x.petsc_vec)
p_.x.scatter_forward()
div_u_s = Function(Q)
div_u_s.interpolate(
Expression(
div(u_s),
Q.element.interpolation_points()
)
)
p_.x.petsc_vec.axpy(-mu, div_u_s.x.petsc_vec)
p_.x.scatter_forward()
The all code is
import gmsh
import meshio
import numpy as np
import tqdm.autonotebook
from mpi4py import MPI
from petsc4py import PETSc
from basix.ufl import element
from dolfinx.fem import (Constant, Function, functionspace,
assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc, Expression)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
create_vector, create_matrix, set_bc)
from dolfinx.io import (VTXWriter, distribute_entity_data, gmshio)
from ufl import (FacetNormal, Identity, Measure, TestFunction, TrialFunction,
as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym, system)
gmsh.initialize()
L = 2.2
H = 0.41
c_x = c_y = 0.2
r = 0.05
gdim = 2
mesh_comm = MPI.COMM_WORLD
model_rank = 0
if mesh_comm.rank == model_rank:
rectangle = gmsh.model.occ.addRectangle(0, 0, 0, L, H, tag=1)
obstacle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)
if mesh_comm.rank == model_rank:
fluid = gmsh.model.occ.cut([(gdim, rectangle)], [(gdim, obstacle)])
gmsh.model.occ.synchronize()
fluid_marker = 1
if mesh_comm.rank == model_rank:
volumes = gmsh.model.getEntities(dim=gdim)
assert (len(volumes) == 1)
gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker)
gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid")
inlet_marker, outlet_marker, wall_marker, obstacle_marker, top_marker = 2, 3, 4, 5, 6
inflow, outflow, walls, obstacle, top = [], [], [], [], []
if mesh_comm.rank == model_rank:
boundaries = gmsh.model.getBoundary(volumes, oriented=False)
for boundary in boundaries:
center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
if np.allclose(center_of_mass, [0, H / 2, 0]):
inflow.append(boundary[1])
elif np.allclose(center_of_mass, [L, H / 2, 0]):
outflow.append(boundary[1])
elif np.allclose(center_of_mass, [L / 2, 0, 0]):
walls.append(boundary[1])
elif np.allclose(center_of_mass, [L / 2, H, 0]):
top.append(boundary[1])
else:
obstacle.append(boundary[1])
gmsh.model.addPhysicalGroup(1, walls, wall_marker)
gmsh.model.setPhysicalName(1, wall_marker, "Walls")
gmsh.model.addPhysicalGroup(1, inflow, inlet_marker)
gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
gmsh.model.addPhysicalGroup(1, outflow, outlet_marker)
gmsh.model.setPhysicalName(1, outlet_marker, "Outlet")
gmsh.model.addPhysicalGroup(1, obstacle, obstacle_marker)
gmsh.model.setPhysicalName(1, obstacle_marker, "Obstacle")
gmsh.model.addPhysicalGroup(1, top, top_marker)
gmsh.model.setPhysicalName(1, top_marker, "Top")
res_min = r / 3
if mesh_comm.rank == model_rank:
distance_field = gmsh.model.mesh.field.add("Distance")
gmsh.model.mesh.field.setNumbers(distance_field, "EdgesList", obstacle)
threshold_field = gmsh.model.mesh.field.add("Threshold")
gmsh.model.mesh.field.setNumber(threshold_field, "IField", distance_field)
gmsh.model.mesh.field.setNumber(threshold_field, "LcMin", res_min)
gmsh.model.mesh.field.setNumber(threshold_field, "LcMax", 0.25 * H)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMin", r)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMax", 2 * H)
min_field = gmsh.model.mesh.field.add("Min")
gmsh.model.mesh.field.setNumbers(min_field, "FieldsList", [threshold_field])
gmsh.model.mesh.field.setAsBackgroundMesh(min_field)
if mesh_comm.rank == model_rank:
gmsh.option.setNumber("Mesh.Algorithm", 8)
gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
gmsh.option.setNumber("Mesh.RecombineAll", 1)
gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.setOrder(2)
gmsh.model.mesh.optimize("Netgen")
mesh, _, ft = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
ft.name = "Facet markers"
t = 0
T = 1 # Final time
dt = 1 / 1600 # Time step size
num_steps = int(T / dt)
k = Constant(mesh, PETSc.ScalarType(dt))
mu = Constant(mesh, PETSc.ScalarType(0.001)) # Dynamic viscosity
rho = Constant(mesh, PETSc.ScalarType(1)) # Density
v_cg2 = element("Lagrange", mesh.topology.cell_name(), 2, shape=(mesh.geometry.dim, ))
s_cg1 = element("Lagrange", mesh.topology.cell_name(), 1)
V = functionspace(mesh, v_cg2)
Q = functionspace(mesh, s_cg1)
W = functionspace(mesh, s_cg1)
fdim = mesh.topology.dim - 1
# Define boundary conditions
"""
class InletVelocity:
def __init__(self, t):
self.t = t
def __call__(self, x):
values = np.zeros((gdim, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = 4 * 1.5 * np.sin(self.t * np.pi / 8) * x[1] * (0.41 - x[1]) / (0.41**2)
return values
"""
class InletVelocity:
def __call__(self, x):
values = np.zeros((gdim, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = 1.0 * (x[1] / 1.0) ** 0.16
return values
#"""
# Inlet
u_inlet = Function(V)
inlet_velocity = InletVelocity()
u_inlet.interpolate(inlet_velocity)
bcu_inflow = dirichletbc(u_inlet, locate_dofs_topological(V, fdim, ft.find(inlet_marker)))
# Walls
u_nonslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
bcu_walls = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, ft.find(wall_marker)), V)
V1, _ = V.sub(1).collapse()
u_slip = Function(V1)
u_slip.x.array[:] = 0.0
bcu_top = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, ft.find(top_marker)), V)
#bcu_top = dirichletbc(u_slip, locate_dofs_topological((V.sub(1), V1), fdim, ft.find(top_marker)), V.sub(1))
# Obstacle
bcu_obstacle = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, ft.find(obstacle_marker)), V)
bcu = [bcu_inflow, bcu_obstacle, bcu_walls]
# Outlet
bcp_top = dirichletbc(PETSc.ScalarType(0), locate_dofs_topological(Q, fdim, ft.find(top_marker)), Q)
bcp_outlet = dirichletbc(PETSc.ScalarType(0), locate_dofs_topological(Q, fdim, ft.find(outlet_marker)), Q)
bcp = [bcp_top, bcp_outlet]
u = TrialFunction(V)
v = TestFunction(V)
u_ = Function(V)
u_.name = "u"
u_s = Function(V)
u_n = Function(V)
u_n1 = Function(V)
p = TrialFunction(Q)
q = TestFunction(Q)
p_ = Function(Q)
p_.name = "p"
phi = Function(Q)
k_n = Function(W); k_n.name = "k"
e_n = Function(W); e_n.name = "e"
u_ref = 2.0
C_mu = 0.09
l_0 = 0.1
k_0 = 0.01 * u_ref**2
e_0 = C_mu**0.75 * (k_0**1.5) / l_0
k_n.interpolate(lambda x: np.full(x.shape[1], k_0))
e_n.interpolate(lambda x: np.full(x.shape[1], e_0))
k_n.x.scatter_forward()
e_n.x.scatter_forward()
mu_t = rho * C_mu * k_n**2 / e_n
f = Constant(mesh, PETSc.ScalarType((0, 0)))
F1 = rho / (2.0 * k) * dot(3.0*u -4.0*u_n + u_n1, v) * dx
F1 += inner(dot(2.0 * u_n - u_n1, nabla_grad(u)), v) * dx
F1 += inner(mu * grad(u), grad(v)) * dx - dot(p_, div(v)) * dx
F1 += dot(f, v) * dx
a1 = form(lhs(F1))
L1 = form(rhs(F1))
A1 = create_matrix(a1)
b1 = create_vector(L1)
a2 = form(dot(grad(p), grad(q)) * dx)
L2 = form(-rho * 3.0 / (2.0 * k) * dot(div(u_s), q) * dx)
A2 = assemble_matrix(a2, bcs=bcp)
A2.assemble()
b2 = create_vector(L2)
a3 = form(rho * dot(u, v) * dx)
L3 = form(rho * dot(u_s, v) * dx - (2.0 * k/3.0) * dot(nabla_grad(phi), v) * dx)
A3 = assemble_matrix(a3)
A3.assemble()
b3 = create_vector(L3)
# Solver for step 1
solver1 = PETSc.KSP().create(mesh.comm)
solver1.setOperators(A1)
solver1.setType(PETSc.KSP.Type.BCGS)
pc1 = solver1.getPC()
pc1.setType(PETSc.PC.Type.JACOBI)
# Solver for step 2
solver2 = PETSc.KSP().create(mesh.comm)
solver2.setOperators(A2)
solver2.setType(PETSc.KSP.Type.MINRES)
pc2 = solver2.getPC()
pc2.setType(PETSc.PC.Type.HYPRE)
pc2.setHYPREType("boomeramg")
# Solver for step 3
solver3 = PETSc.KSP().create(mesh.comm)
solver3.setOperators(A3)
solver3.setType(PETSc.KSP.Type.CG)
pc3 = solver3.getPC()
pc3.setType(PETSc.PC.Type.SOR)
from pathlib import Path
folder = Path("/home/gerard/Desktop/benchmark_rotational")
folder.mkdir(exist_ok=True, parents=True)
vtx_u = VTXWriter(mesh.comm, folder / "dfg2D-3-u.bp", [u_], engine="BP4")
vtx_p = VTXWriter(mesh.comm, folder / "dfg2D-3-p.bp", [p_], engine="BP4")
vtx_u.write(t)
vtx_p.write(t)
progress = tqdm.autonotebook.tqdm(desc="Solving PDE", total=num_steps)
for i in range(num_steps):
progress.update(1)
# Update current time step
t += dt
# Update inlet velocity
inlet_velocity.t = t
u_inlet.interpolate(inlet_velocity)
# Step 1: Tentative velocity step
A1.zeroEntries()
assemble_matrix(A1, a1, bcs=bcu)
A1.assemble()
with b1.localForm() as loc:
loc.set(0)
assemble_vector(b1, L1)
apply_lifting(b1, [a1], [bcu])
b1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b1, bcu)
solver1.solve(b1, u_s.x.petsc_vec)
u_s.x.scatter_forward()
# Step 2: Pressure corrrection step
with b2.localForm() as loc:
loc.set(0)
assemble_vector(b2, L2)
apply_lifting(b2, [a2], [bcp])
b2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b2, bcp)
solver2.solve(b2, phi.x.petsc_vec)
phi.x.scatter_forward()
p_.x.petsc_vec.axpy(1, phi.x.petsc_vec)
p_.x.scatter_forward()
div_u_s = Function(Q)
div_u_s.interpolate(
Expression(
div(u_s),
Q.element.interpolation_points()
)
)
p_.x.petsc_vec.axpy(-mu, div_u_s.x.petsc_vec)
p_.x.scatter_forward()
# Step 3: Velocity correction step
with b3.localForm() as loc:
loc.set(0)
assemble_vector(b3, L3)
b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
solver3.solve(b3, u_.x.petsc_vec)
u_.x.scatter_forward()
# Write solutions to file
vtx_u.write(t)
vtx_p.write(t)
# Update variable with solution form this time step
with u_.x.petsc_vec.localForm() as loc_, u_n.x.petsc_vec.localForm() as loc_n, u_n1.x.petsc_vec.localForm() as loc_n1:
loc_n.copy(loc_n1)
loc_.copy(loc_n)
progress.close()
vtx_u.close()
vtx_p.close()