My intension is to calculate the shear stresses and normal stresses on a pipe under Poiseuille flow (Test problem 1: Channel flow )
I found a FEniCS code that have done the same thing, but now I’m having some trouble to convert it into the FEniCSx syntax.
The FEniCS code is as follows: (Source here)
And I was able to convert upto the following:
Required imports (including the items necessary for the test problem):
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import pyvista
from dolfinx.fem import Constant,Function, FunctionSpace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_unit_square
from dolfinx.plot import create_vtk_mesh
from ufl import (FacetNormal,FacetArea,as_vector, FiniteElement, Identity,TestFunction, TrialFunction, VectorElement,
div, dot, ds, dx, inner, lhs, nabla_grad, rhs, sym,inner)
My attempt to convert the FEniCS code into the FEniCSx syntax:
def epsilon(u):
return sym(nabla_grad(u))
def sigma(u, p):
return 2*mu*epsilon(u) - p*Identity(u.geometric_dimension())
#The above two values will anyway be calculated in the test problem.
n = FacetNormal(mesh)
T = -sigma(u, p)*n
Tn = inner(T, n) # Inner product: https://jorgensd.github.io/dolfinx-tutorial/chapter1/fundamentals_code.html?highlight=inner%20product#expressing-inner-products
Tt = T - Tn*n
#I'm changing the piecewise constant to, order 1 and order 2
p_cg1 = FiniteElement('CG', mesh.ufl_cell(), 1)
v_cg2 = VectorElement('CG', mesh.ufl_cell(), 2)
scalar = FunctionSpace(mesh, p_cg1)
vector = FunctionSpace(mesh, v_cg2)
v = TestFunction(scalar)
w = TestFunction(vector)
normal_stress = Function(scalar)
shear_stress = Function(vector)
Ln = (1/FacetArea(mesh))* v*Tn * ds
Lt = (1/FacetArea(mesh))* inner(w,Tt) * ds
But, then I’m confused on how to properly assemble it and obtain the shear stresses on the boundary (We have a unit square as the domain in the Test problem 1)
Can you please help me to proceed after this step.
Thank you in advance
A minimum reproducible example for the NS problem is:
(The exact code in the Test problem 1: Channel flow )
from mpi4py import MPI from petsc4py import PETSc import numpy as np import pyvista from dolfinx.fem import Constant,Function, FunctionSpace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc from dolfinx.io import XDMFFile from dolfinx.mesh import create_unit_square from dolfinx.plot import create_vtk_mesh from ufl import (FacetNormal,FacetArea,as_vector, FiniteElement, Identity,TestFunction, TrialFunction, VectorElement, div, dot, ds, dx, inner, lhs, nabla_grad, rhs, sym,inner) mesh = create_unit_square(MPI.COMM_WORLD, 10, 10) t = 0 T = 10 num_steps = 500 dt = T/num_steps v_cg2 = VectorElement("CG", mesh.ufl_cell(), 2) s_cg1 = FiniteElement("CG", mesh.ufl_cell(), 1) V = FunctionSpace(mesh, v_cg2) Q = FunctionSpace(mesh, s_cg1) u = TrialFunction(V) v = TestFunction(V) p = TrialFunction(Q) q = TestFunction(Q) def walls(x): return np.logical_or(np.isclose(x[1],0), np.isclose(x[1],1)) wall_dofs = locate_dofs_geometrical(V, walls) u_noslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType) bc_noslip = dirichletbc(u_noslip, wall_dofs, V) def inflow(x): return np.isclose(x[0], 0) inflow_dofs = locate_dofs_geometrical(Q, inflow) bc_inflow = dirichletbc(PETSc.ScalarType(8), inflow_dofs, Q) def outflow(x): return np.isclose(x[0], 1) outflow_dofs = locate_dofs_geometrical(Q, outflow) bc_outflow = dirichletbc(PETSc.ScalarType(0), outflow_dofs, Q) bcu = [bc_noslip] bcp = [bc_inflow, bc_outflow] u_n = Function(V) u_n.name = "u_n" U = 0.5 * (u_n + u) n = FacetNormal(mesh) f = Constant(mesh, PETSc.ScalarType((0,0))) k = Constant(mesh, PETSc.ScalarType(dt)) mu = Constant(mesh, PETSc.ScalarType(1)) rho = Constant(mesh, PETSc.ScalarType(1)) # Define strain-rate tensor def epsilon(u): return sym(nabla_grad(u)) # Define stress tensor def sigma(u, p): return 2*mu*epsilon(u) - p*Identity(u.geometric_dimension()) # https://jorgensd.github.io/dolfinx-tutorial/chapter2/navierstokes.html # Define the variational problem for the first step p_n = Function(Q) p_n.name = "p_n" F1 = rho*dot((u - u_n) / k, v)*dx F1 += rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx F1 += inner(sigma(U, p_n), epsilon(v))*dx F1 += dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds F1 -= dot(f, v)*dx a1 = form(lhs(F1)) L1 = form(rhs(F1)) A1 = assemble_matrix(a1, bcs=bcu) A1.assemble() b1 = create_vector(L1) # Define variational problem for step 2 u_ = Function(V) a2 = form(dot(nabla_grad(p), nabla_grad(q))*dx) L2 = form(dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx) A2 = assemble_matrix(a2, bcs=bcp) A2.assemble() b2 = create_vector(L2) # Define variational problem for step 3 p_ = Function(Q) a3 = form(dot(u, v)*dx) L3 = form(dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), 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.HYPRE) pc1.setHYPREType("boomeramg") # Solver for step 2 solver2 = PETSc.KSP().create(mesh.comm) solver2.setOperators(A2) solver2.setType(PETSc.KSP.Type.BCGS) 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) xdmf = XDMFFile(mesh.comm, "poiseuille.xdmf", "w") xdmf.write_mesh(mesh) xdmf.write_function(u_n, t) xdmf.write_function(p_n, t) import numpy as np def u_exact(x): values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType) values[0] = 4*x[1]*(1.0 - x[1]) return values u_ex = Function(V) u_ex.interpolate(u_exact) L2_error = form(dot(u_ - u_ex, u_ - u_ex)*dx) for i in range(num_steps): # Update current time step t += dt # Step 1: Tentative veolcity step with b1.localForm() as loc_1: loc_1.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_.vector) u_.x.scatter_forward() # Step 2: Pressure corrrection step with b2.localForm() as loc_2: loc_2.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, p_.vector) p_.x.scatter_forward() # Step 3: Velocity correction step with b3.localForm() as loc_3: loc_3.set(0) assemble_vector(b3, L3) b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) solver3.solve(b3, u_.vector) u_.x.scatter_forward() # Update variable with solution form this time step u_n.x.array[:] = u_.x.array[:] p_n.x.array[:] = p_.x.array[:] # Write solutions to file xdmf.write_function(u_n, t) xdmf.write_function(p_n, t) # Compute error at current time-step error_L2 = np.sqrt(mesh.comm.allreduce(assemble_scalar(L2_error), op=MPI.SUM)) error_max = mesh.comm.allreduce(np.max(u_.vector.array - u_ex.vector.array), op=MPI.MAX) # Print error only every 20th step and at the last step # if (i % 20 == 0) or (i == num_steps - 1): # print(f"Time {t:.2f}, L2-error {error_L2:.2e}, Max error {error_max:.2e}") # Close xmdf file xdmf.close()To plot the velocity field and to make sure it works fine:
# pyvista.set_jupyter_backend("pythreejs") # pyvista.set_jupyter_backend("ipygany") pyvista.set_jupyter_backend("ipyvtklink") pyvista.start_xvfb() topology, cell_types, geometry = create_vtk_mesh(V) values = np.zeros((geometry.shape[0], 3), dtype=np.float64) values[:, :len(u_n)] = u_n.x.array.real.reshape((geometry.shape[0], len(u_n))) # Create a point cloud of glyphs function_grid = pyvista.UnstructuredGrid(topology, cell_types, geometry) function_grid["u"] = values glyphs = function_grid.glyph(orient="u", factor=0.2) # Create a pyvista-grid for the mesh grid = pyvista.UnstructuredGrid(*create_vtk_mesh(mesh, mesh.topology.dim)) # Create plotter plotter = pyvista.Plotter() # plotter.add_mesh(grid, style="wireframe", color="k") plotter.add_mesh(glyphs) plotter.view_xy() if not pyvista.OFF_SCREEN: plotter.show() else: pyvista.start_xvfb() fig_as_array = plotter.screenshot("glyphs.png")The code that I have so far, for the shear stress calculation:
n = FacetNormal(mesh) T = -sigma(u, p)*n Tn = inner(T, n) # Inner product: https://jorgensd.github.io/dolfinx-tutorial/chapter1/fundamentals_code.html?highlight=inner%20product#expressing-inner-products Tt = T - Tn*n #I'm changing the piecewise constant to, order 1 and order 2 p_cg1 = FiniteElement('CG', mesh.ufl_cell(), 1) v_cg2 = VectorElement('CG', mesh.ufl_cell(), 2) scalar = FunctionSpace(mesh, p_cg1) vector = FunctionSpace(mesh, v_cg2) v = TestFunction(scalar) w = TestFunction(vector) normal_stress = Function(scalar) shear_stress = Function(vector) Ln = (1/FacetArea(mesh))* v*Tn * ds Lt = (1/FacetArea(mesh))* inner(w,Tt) * ds assemble_vector(normal_stress.vector, form(Ln)) assemble_vector(normal_stress.vector, form(Lt))