# %% # Advection–Diffusion or Fluid Transport in a 2D Periodic Channel # =============================================================== # # This example illustrates how to solve a 2D periodic channel problem # using either the advection–diffusion equation (ADV_DIFF = 1) # or the incompressible Navier–Stokes equations (ADV_DIFF = 0). # # The domain is a rectangular channel with periodic lateral boundaries # and Dirichlet conditions at the top and bottom. # # Keywords # -------- # Advection-Diffusion, Periodic Boundaries, Mesh Adaptation import os, shutil, gmsh, sys import numpy as np from migflow import fluid, advdiff, gmsh_io # %% # Output Directory # ---------------- # Create a clean output directory for simulation results. outputdir = "output_2d_marker_advection" if len(sys.argv) < 2 else sys.argv[1] shutil.rmtree(outputdir, ignore_errors=True) os.makedirs(outputdir) # %% # Problem Configuration # --------------------- ADV_DIFF = 1 # 1 → advection–diffusion, 0 → fluid simulation # %% # Helper Class for Mesh Sizing # ---------------------------- class MeshSize: """Callable mesh-size function for gmsh.""" def __init__(self, lc): self._lc = lc def __call__(self, dim, tag, x, y, z, lc): return self._lc # %% # Geometrical parameters # ---------------------- width, height = 4.0, 1.0 origin = [-width / 2, -height / 2] lc = 0.1 mesh_size = MeshSize(lc) # %% # Physical Parameters # ------------------- g = np.array([0, 0]) # gravity (ignored for adv-diff) rho = 1.204 # nitrogen density [kg/m³] nu = 1.8e-5 # kinematic viscosity [m²/s] mu = rho * nu # dynamic viscosity [kg/(m·s)] alpha = 3e-3 # thermal dilatation coefficient k = 2e-2 # thermal conductivity [W/mK] cp = 1010e-3 # heat capacity [kJ/kg·K] # %% # Generate Initial Mesh # --------------------- def gen_mesh(width, height, mesh_size, origin=np.array([0, 0])): """Generate a rectangular periodic mesh in gmsh. Parameters ---------- width : float Domain width. height : float Domain height. mesh_size : MeshSize Callable object controlling element size. origin : np.ndarray, optional Origin of the rectangle. """ origin = np.asarray(origin) gmsh.model.add("box") gmsh.model.occ.add_rectangle(origin[0], origin[1], 0, width, height) gmsh.model.occ.synchronize() def get_line(x0, x1, eps=1e-6): r = gmsh.model.get_entities_in_bounding_box( x0[0] - eps, x0[1] - eps, -eps, x1[0] + eps, x1[1] + eps, eps, 1 ) return [tag for _, tag in r] h, w = height, width l_tag = get_line(origin + [0, 0], origin + [0, h]) r_tag = get_line(origin + [w, 0], origin + [w, h]) # Define periodicity (right ↔ left) affine_matrix = np.eye(4, 4, dtype=np.float64) affine_matrix[0, 3] = width gmsh.model.mesh.set_periodic(1, r_tag, l_tag, affine_matrix.reshape(-1)) gmsh.model.occ.synchronize() # Define physical groups gmsh.model.add_physical_group( 1, get_line(origin + [0, 0], origin + [w, 0]), name="Bottom" ) gmsh.model.add_physical_group( 1, get_line(origin + [0, h], origin + [w, h]), name="Top" ) gmsh.model.add_physical_group(2, [1], name="domain") gmsh.model.mesh.set_size_callback(mesh_size) gmsh.initialize() gen_mesh(width, height, mesh_size, origin) gmsh.model.mesh.generate(2) # %% # Problem Definition # ------------------ def define_problem(adv_diff): """Create either an advection–diffusion or fluid problem.""" if adv_diff: # --- Advection–Diffusion problem --- # Discontinuous P1, the scheme this case was written for: the kernel # now carries the interior-face upwind flux that couples one element to # the next, so a cg=False field transports (it could not before), and # the slope limiter keeps it monotone. Measured here: the marker # translates exactly at u t, its integral is conserved to machine # precision, and it stays in [0, 1]. c = advdiff.AdvDiffProblem2(0, 0, 1.0, explicit_scheme=True, cg=False) gmsh_io.load_msh(c, None) c.set_weak_boundary("Bottom") c.set_weak_boundary("Top") velocity = c.velocity().get() velocity[:, 0] = 1.0 # uniform flow to the right c.velocity().set_from_host(velocity) # The solution is discontinuous, so its dofs are per element and per # local node; scatter the node coordinates onto them. s0 = c.solution().get() x = np.empty((s0.shape[0], 3)) x[c.get_mapping(c.get_solution_element())] = c.coordinates()[c.elements()] s0[x[:, 0] > -1e-8] = 1 # initialize concentration at inflow c.solution().set_from_host(s0) return c else: # --- Fluid problem (reference) --- f = fluid.FluidProblem2(g, mu, rho, density_element="triangle_p1dg", viscosity_element="triangle_p1dg") gmsh_io.load_msh(f, None) f.set_mean_pressure(0) f.set_wall_boundary("Bottom") f.set_wall_boundary("Top") s0 = np.zeros_like(f.solution().get()) s0[f.velocity_index()[:, 0]] = 1.0 # initial velocity in x-direction f.solution().set_from_host(s0) c0 = np.zeros_like(f.coordinates()[:, 0]) xc0 = f.coordinates() c0[xc0[:, 0] > -1e-8] = 1 f.set_concentration_cg(c0) return f p = define_problem(ADV_DIFF) # %% # Mesh Adaptation # --------------- mesh_size = MeshSize(0.05) gen_mesh(width, height, mesh_size, origin) gmsh.model.mesh.clear() gmsh.model.mesh.generate(2) if ADV_DIFF: gmsh_io.adapt_msh(p, None) velocity = p.velocity().get() velocity[:, 0] = 1.0 p.velocity().set_from_host(velocity) else: gmsh_io.adapt_msh(p, None) # %% # Simulation Parameters # --------------------- dt = 0.01 # time step [s] tEnd = 1.0 # final time [s] outf = 1 # output frequency # %% # Time Integration Loop # --------------------- t, i = 0, 0 while t < tEnd: print(f"{i:4d}, {t:.6g}/{tEnd:.6g}") if i % outf == 0: p.write_mig(outputdir, t) if ADV_DIFF: p.explicit_euler(dt, limiter=True) else: p.advance_concentration(dt) t += dt i += 1 # %% # Plot # ---- # .. code-block:: shell # # python3 -m migflow.plot.migplot output_2d_marker_advection --actors fluid --problem-type advection-diffusion --fluid-field solution --fluid-vmin 0 --fluid-vmax 1 --show-edges 1 --element-type triangle_p1dg # %% # Artifacts # --------- # - output_2d_marker_advection/animation.mp4