Hi,
I have this attached geometry consisted of one inlet and two outlets, but it keeps diverging whatever I tried to do (enhancing the mesh resolution, decreasing the time step, decreasing the inlet flow rate, using different preconditioner and stabilization) one thing I noticed as well is that I’m applying the parabolic inlet profile but the region after the inlet suddenly doesn’t follow that (kinda discontnuity) and after a few steps the velocity explode and backflow happen.
Geometry:
The inlet profile discontinuity:
My script
#!/usr/bin/env python3
# coding: utf-8
“”"
Optimized IPCS (Oasis/DFG-style) scheme in classic FEniCS (dolfin):
- BDF2 in time: (3 u^{n+1} - 4 u^n + u^{n-1}) / (2 dt) - Convection in A1 with explicit convecting velocity: u_conv = 2 u^n - u^{n-1} i.e. (u_conv · ∇ u^{n+1}, v) - Diffusion implicit: mu (∇u^{n+1}, ∇v) - Pressure term in step 1 uses current pressure p^n (held in p_) - Pressure correction Poisson: (∇phi, ∇q) = -(rho * 3/(2 dt)) (div u_s, q) - Rotational pressure update: p^{n+1} = p^n + phi - mu * div(u_s) (div(u_s) is projected into Q via scalar mass solve) - Velocity correction: rho (u^{n+1}, v) = rho (u_s, v) - (2 dt/3) (∇phi, v) Performance-minded: - Preassemble constant matrices (vector mass Mv, stiffness K, scalar mass Ms, A2, A3) - Reuse convection matrix storage (assemble into tensor=C each step) - Reuse RHS vectors and temporary vectors - Print AND save every single step (as you requested) Mesh input: Mesh/Volume_mesh.xml Mesh/BCnodeFacets.xml Boundary tags: 1 inlet, 2 wall, 3/4 outlets """ from __future__ import print_function import os import time import numpy as np from fenics import * # ============================================================================== # MPI / LOGGING # ============================================================================== comm = MPI.comm_world rank = MPI.rank(comm) def log(msg): if rank == 0: print(msg, flush=True) def banner(msg): if rank == 0: print("\n" + "=" * 80) print(msg) print("=" * 80 + "\n", flush=True) # ============================================================================== # FEniCS SETTINGS # ============================================================================== parameters["ghost_mode"] = "shared_facet" parameters["form_compiler"]["optimize"] = True parameters["form_compiler"]["quadrature_degree"] = 2 set_log_level(LogLevel.PROGRESS) # ============================================================================== # USER PARAMETERS # ============================================================================== banner("USER PARAMETERS") MU = 0.04 RHO = 1.06 INLET_ID = 1 WALL_ID = 2 OUTLET_IDS = [3, 4] DT = 1e-5 T_MAX = 1.0 MAX_STEPS = int(np.ceil(T_MAX / DT)) INLET_U_TARGET = 57.0 N_RAMP = 1000 OUT_DIR = "results_bdf2_rotational_ipcs" os.makedirs(OUT_DIR, exist_ok=True) U_XDMF = os.path.join(OUT_DIR, "velocity.xdmf") P_XDMF = os.path.join(OUT_DIR, "pressure.xdmf") log(f"MU={MU}, RHO={RHO}") log(f"DT={DT}, T_MAX={T_MAX}, MAX_STEPS={MAX_STEPS}") log(f"INLET_U_TARGET={INLET_U_TARGET}, N_RAMP={N_RAMP}") log(f"OUT_DIR={OUT_DIR}") # ============================================================================== # LOAD MESH + FACET TAGS # ============================================================================== banner("LOADING MESH") mesh = Mesh("Mesh_1k/Volume_mesh.xml") facets = MeshFunction("size_t", mesh, "Mesh_1k/BCnodeFacets.xml") if rank == 0: log(f"[MESH] Cells : {mesh.num_cells()}") log(f"[MESH] Vertices : {mesh.num_vertices()}") log(f"[MESH] Dim : {mesh.geometry().dim()}") log(f"[FACETS] Unique IDs: {np.unique(facets.array())}") # ============================================================================== # FUNCTION SPACES # ============================================================================== banner("FUNCTION SPACES") V = VectorFunctionSpace(mesh, "CG", 1) Q = FunctionSpace(mesh, "CG", 1) u = TrialFunction(V) v = TestFunction(V) p = TrialFunction(Q) q = TestFunction(Q) # velocity states u_nm1 = Function(V, name="u_nm1") # u^{n-1} u_n = Function(V, name="u_n") # u^{n} u_s = Function(V, name="u_star") # tentative u_new = Function(V, name="u") # u^{n+1} # pressure states p_ = Function(Q, name="p_n") # p^{n} p_new = Function(Q, name="p") # p^{n+1} phi_ = Function(Q, name="phi") # pressure correction # convecting velocity for step 1 u_conv = Function(V, name="u_conv") log(f"[DoF] V.dim = {V.dim()} | Q.dim = {Q.dim()}") # ============================================================================== # INLET PROFILE (your PCA normal version from earlier, kept compatible) # ============================================================================== banner("INLET PROFILE") # If you want your PCA-normal parabolic inlet, keep it. # Here we reuse your earlier idea: infer inlet patch vertices and compute PCA normal + radius. inlet_facet_ids = np.where(facets.array() == INLET_ID)[0] if inlet_facet_ids.size == 0: raise RuntimeError("No inlet facets found. Check INLET_ID and facet tags.") vertex_ids = [] for fid in inlet_facet_ids: vertex_ids.extend(Facet(mesh, int(fid)).entities(0)) vertex_ids = np.unique(np.array(vertex_ids, dtype=np.int64)) coords = mesh.coordinates()[vertex_ids] inlet_center = coords.mean(axis=0) X = coords - inlet_center Ccov = np.dot(X.T, X) / max(X.shape[0], 1) eigvals, eigvecs = np.linalg.eigh(Ccov) normal = eigvecs[:, np.argmin(eigvals)] normal = normal / np.linalg.norm(normal) domain_center = mesh.coordinates().mean(axis=0) if np.dot(normal, domain_center - inlet_center) < 0: normal *= -1 radius = float(np.max(np.linalg.norm(coords - inlet_center, axis=1))) if rank == 0: log(f"[INLET] center = {inlet_center}") log(f"[INLET] radius = {radius}") log(f"[INLET] normal = {normal}") class ParabolicInlet(UserExpression): def __init__(self, center, radius, normal, Umax, **kwargs): self.c = np.array(center, dtype=np.float64) self.R = float(radius) self.n = np.array(normal, dtype=np.float64) self.Umax = float(Umax) super().__init__(**kwargs) def value_shape(self): return (mesh.geometry().dim(),) def eval(self, values, x): d = np.array(x, dtype=np.float64) - self.c rvec = d - np.dot(d, self.n) * self.n r = np.linalg.norm(rvec) if r <= self.R: values[:] = self.Umax * (1.0 - (r / self.R) ** 2) * self.n else: values[:] = 0.0 inlet_expr = ParabolicInlet( center=inlet_center, radius=radius, normal=normal, Umax=0.0, degree=2 ) # ============================================================================== # BOUNDARY CONDITIONS # ============================================================================== banner("BOUNDARY CONDITIONS") bcu = [ DirichletBC(V, inlet_expr, facets, INLET_ID), DirichletBC(V, Constant((0.0, 0.0, 0.0)), facets, WALL_ID), ] # Pressure Dirichlet on outlets (same style as your earlier scripts)press bcp = [DirichletBC(Q, Constant(0.0), facets, i) for i in OUTLET_IDS] log(f"[BC] inlet={INLET_ID}, wall={WALL_ID}, outlets={OUTLET_IDS}") # ============================================================================== # FORMS (BDF2 + explicit convecting vel, rotational correction) # ============================================================================== banner("FORMS") dt = Constant(DT) rho = Constant(RHO) mu = Constant(MU) nrm = FacetNormal(mesh) # Constant matrices / patterns mass_v_form = inner(u, v) * dx stiff_v_form = inner(grad(u), grad(v)) * dx mass_s_form = inner(p, q) * dx poisson_form = inner(grad(p), grad(q)) * dx mass_rho_form = rho * inner(u, v) * dx # Convection bilinear form with explicit convecting velocity u_conv: # (u_conv · ∇ u^{n+1}, v) conv_form = inner(dot(u_conv, nabla_grad(u)), v) * dx # Pressure term in step 1 uses p^n: # - (p^n, div(v)) dx + <p^n n, v> ds ds = Measure("ds", domain=mesh, subdomain_data=facets) press_rhs_form = (- p_ * div(v) * dx) + dot(p_ * nrm, v) * (ds(3) + ds(4)) # press_rhs_form = (- p_ * div(v) * dx) + dot(p_ * nrm, v) * ds # Step 2 RHS: # (∇phi, ∇q) = -(rho * 3/(2dt)) (div u_s, q) rhs_phi_form = -(rho * (3.0 / (2.0 * dt))) * div(u_s) * q * dx # Velocity correction RHS: # rho (u^{n+1}, v) = rho (u_s, v) - (2dt/3) (∇phi, v) velcorr_rhs_form = rho * inner(u_s, v) * dx - (2.0 * dt / 3.0) * inner(grad(phi_), v) * dx # For rotational pressure update: project div(u_s) into Q: # Solve: (d, q) = (div(u_s), q) div_rhs_form = div(u_s) * q * dx # ============================================================================== # PREASSEMBLE CONSTANT MATRICES # ============================================================================== banner("ASSEMBLING CONSTANT MATRICES") Mv = assemble(mass_v_form) # vector mass K = assemble(stiff_v_form) # vector stiffness Ms = assemble(mass_s_form) # scalar mass A2 = assemble(poisson_form) # Poisson matrix A3 = assemble(mass_rho_form) # rho * mass # Apply pressure Dirichlet to Poisson matrix once for bc in bcp: bc.apply(A2) # Assemble and lock Mv.ident_zeros() K.ident_zeros() Ms.ident_zeros() A2.ident_zeros() A3.ident_zeros() # Convection matrix placeholder (same sparsity as Mv typically for CG1 velocity) C = PETScMatrix() assemble(conv_form, tensor=C) # allocate pattern once C.zero() # Tentative velocity matrix placeholder A1 = PETScMatrix() assemble(mass_v_form, tensor=A1) # allocate once A1.zero() # Reusable vectors b1 = Vector(); Mv.init_vector(b1, 0) tmp = Vector(); Mv.init_vector(tmp, 0) tmp2 = Vector(); Mv.init_vector(tmp2, 0) b2 = Vector(); A2.init_vector(b2, 0) b3 = Vector(); A3.init_vector(b3, 0) bdiv = Vector(); Ms.init_vector(bdiv, 0) # ============================================================================== # SOLVERS # ============================================================================== banner("SOLVERS") # Step 1: tentative velocity solver_u = KrylovSolver("gmres", "ilu") solver_u.parameters["relative_tolerance"] = 1e-8 solver_u.parameters["absolute_tolerance"] = 1e-10 solver_u.parameters["maximum_iterations"] = 1000 solver_u.parameters["monitor_convergence"] = False # Step 2: pressure correction solver_phi = KrylovSolver("cg", "hypre_amg") solver_phi.parameters["relative_tolerance"] = 1e-8 solver_phi.parameters["absolute_tolerance"] = 1e-10 solver_phi.parameters["maximum_iterations"] = 2000 solver_phi.parameters["monitor_convergence"] = False # Step 3: velocity correction (rho mass) solver_corr = KrylovSolver("cg", "jacobi") solver_corr.parameters["relative_tolerance"] = 1e-8 solver_corr.parameters["absolute_tolerance"] = 1e-10 solver_corr.parameters["maximum_iterations"] = 500 solver_corr.parameters["monitor_convergence"] = False # Rotational projection: scalar mass solve for div(u_s) solver_div = KrylovSolver("cg", "jacobi") solver_div.parameters["relative_tolerance"] = 1e-8 solver_div.parameters["absolute_tolerance"] = 1e-10 solver_div.parameters["maximum_iterations"] = 500 solver_div.parameters["monitor_convergence"] = False # ============================================================================== # OUTPUT (SAVE EVERY STEP) # ============================================================================== banner("OUTPUT") xdmf_u = XDMFFile(comm, U_XDMF) xdmf_p = XDMFFile(comm, P_XDMF) xdmf_u.parameters["flush_output"] = True xdmf_u.parameters["rewrite_function_mesh"] = False xdmf_p.parameters["flush_output"] = True xdmf_p.parameters["rewrite_function_mesh"] = False # ============================================================================== # INITIALIZATION # ============================================================================== banner("INITIALIZATION") u_nm1.vector().zero() u_n.vector().zero() p_.vector().zero() t = 0.0 start_time = time.time() # Write initial condition u_new.assign(u_n) p_new.assign(p_) xdmf_u.write(u_new, t) xdmf_p.write(p_new, t) # ============================================================================== # TIME LOOP (PRINT & SAVE EVERY STEP) # ============================================================================== banner("TIME LOOP (BDF2 + ROTATIONAL IPCS)") for step in range(MAX_STEPS): # Update time and inlet ramp t = (step + 1) * DT inlet_expr.Umax = INLET_U_TARGET * (step + 1) / N_RAMP if step < N_RAMP else INLET_U_TARGET # Build convecting velocity: # u_conv = 2 u^n - u^{n-1} (BDF2 explicit) u_conv.vector().zero() u_conv.vector().axpy(2.0, u_n.vector()) u_conv.vector().axpy(-1.0, u_nm1.vector()) # Assemble convection matrix into preallocated tensor C.zero() assemble(conv_form, tensor=C) # ----------------------------- # STEP 1: Tentative velocity # # For step==0, do BDF1 startup: # (rho/dt)(u - u^n) + (u^n · ∇u) - div(mu ∇u) + ∇p^n = 0 # For step>=1, BDF2: # (rho/(2dt))(3u - 4u^n + u^{n-1}) + (u_conv · ∇u) - div(mu ∇u) + ∇p^n = 0 # ----------------------------- A1.zero() b1.zero() if step == 0: # Matrix: (rho/dt) M + mu K + C A1.axpy(float(RHO / DT), Mv, True) A1.axpy(float(MU), K, True) A1.axpy(1.0, C, True) # RHS: (rho/dt) M u^n + pressure RHS Mv.mult(u_n.vector(), b1) b1 *= float(RHO / DT) else: # Matrix: (3 rho / (2 dt)) M + mu K + C A1.axpy(float(3.0 * RHO / (2.0 * DT)), Mv, True) A1.axpy(float(MU), K, True) A1.axpy(1.0, C, True) # RHS: (rho/(2dt)) M (4 u^n - u^{n-1}) + pressure RHS # tmp = 4 u^n - u^{n-1} tmp2.zero() tmp2.axpy(4.0, u_n.vector()) tmp2.axpy(-1.0, u_nm1.vector()) Mv.mult(tmp2, b1) b1 *= float(RHO / (2.0 * DT)) # Add pressure RHS (assembled vector) b_press = assemble(press_rhs_form) b1.axpy(1.0, b_press) # Apply velocity Dirichlet BCs for bc in bcu: bc.apply(A1, b1) # Solve for tentative velocity solver_u.solve(A1, u_s.vector(), b1) # ----------------------------- # STEP 2: Pressure correction # A2 phi = -(rho * 3/(2dt)) div(u_s) # ----------------------------- b2.zero() assemble(rhs_phi_form, tensor=b2) # Apply pressure Dirichlet BCs to RHS for bc in bcp: bc.apply(b2) solver_phi.solve(A2, phi_.vector(), b2) # Pressure update: p = p + phi p_new.vector().zero() p_new.vector().axpy(1.0, p_.vector()) p_new.vector().axpy(1.0, phi_.vector()) # Rotational correction: p = p - mu * div(u_s) # Compute d in Q such that (d, q) = (div(u_s), q) bdiv.zero() assemble(div_rhs_form, tensor=bdiv) d = Function(Q, name="div_u_s_proj") for bc in bcp: bc.apply(Ms, bdiv) # safe: keeps pressure dirichlet consistent at outlets if you want solver_div.solve(Ms, d.vector(), bdiv) p_new.vector().axpy(-float(MU), d.vector()) # Enforce outlet pressure strongly (keeps consistent) for bc in bcp: bc.apply(p_new.vector()) # ----------------------------- # STEP 3: Velocity correction # A3 u = rho M u = rho M u_s - (2dt/3) (∇phi, v) # ----------------------------- b3.zero() assemble(velcorr_rhs_form, tensor=b3) # Apply velocity BCs to corrected velocity solve for bc in bcu: bc.apply(A3, b3) solver_corr.solve(A3, u_new.vector(), b3) n = FacetNormal(mesh) ds_in = Measure("ds", domain=mesh, subdomain_data=facets, subdomain_id=INLET_ID) ds_o3 = Measure("ds", domain=mesh, subdomain_data=facets, subdomain_id=3) ds_o4 = Measure("ds", domain=mesh, subdomain_data=facets, subdomain_id=4) Q_in = assemble(dot(u_new, n) * ds_in) Q_o3 = assemble(dot(u_new, n) * ds_o3) Q_o4 = assemble(dot(u_new, n) * ds_o4) if rank == 0: print(f" flux: Qin={Q_in:+.6e} | Qo3={Q_o3:+.6e} | Qo4={Q_o4:+.6e} | sum={Q_in+Q_o3+Q_o4:+.3e}") div_l2 = np.sqrt(assemble((div(u_new)**2) * dx)) if rank == 0: print(f" div_L2={div_l2:.3e}", flush=True) # ----------------------------- # PRINT EVERY STEP # ----------------------------- # relative change measure (against u_n) du = u_new.vector().copy() du.axpy(-1.0, u_n.vector()) rel = du.norm("l2") / max(u_new.vector().norm("l2"), 1e-30) # max velocity magnitude (approx, using l_inf of vector dofs) u_max = u_new.vector().norm("linf") if rank == 0: print( f"[STEP {step:7d}] t={t:.6e} | Uin={inlet_expr.Umax:9.3f} | rel={rel:.3e} | u_max={u_max:.3e}", flush=True ) # ----------------------------- # SAVE EVERY STEP # ----------------------------- xdmf_u.write(u_new, t) xdmf_p.write(p_new, t) # Roll time levels: u^{n-1} <- u^n, u^n <- u^{n+1}, p^n <- p^{n+1} u_nm1.assign(u_n) u_n.assign(u_new) p_.assign(p_new) elapsed = time.time() - start_time banner("DONE") log(f"Total wall time = {elapsed:.2f} s")
and this is my mesh