I’ve been working on this problem for a couple of weeks now but, I couldn’t manage to get a proper answer.
My goal is to compute the Shear stress and the normal stress on a 2D Channel flow.
Lets start with a sample code for the Navier stokes solver on a unit square: It might look long but it is just the exact replica from the FEniCSx tutorial Test problem 1: Channel flow :
import numpy as np
import gmsh
import sys
from dolfinx import plot
# Initialize gmsh:
from mpi4py import MPI
#from dolfinx import mesh
import pyvista
from dolfinx.io.gmshio import model_to_mesh
import matplotlib.pyplot as plt
from dolfinx import geometry
# from dolfinx import fem,mesh
import ufl
from petsc4py import PETSc
from prettytable import PrettyTable
import random
from dolfinx.mesh import create_unit_square,exterior_facet_indices
from dolfinx.cpp.mesh import entities_to_geometry
from ufl import (FacetNormal, FiniteElement, FacetArea,Identity,TestFunction, TrialFunction, VectorElement,
div, dot, ds, dx, inner, lhs, nabla_grad, rhs, sym,as_vector,grad,SpatialCoordinate)
from dolfinx.fem import Expression,Constant,Function, FunctionSpace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical,locate_dofs_topological
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
from dolfinx.io import XDMFFile
from dolfinx.plot import create_vtk_mesh
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())
# 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)
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)
# Close xmdf file
xdmf.close()
Now we have the velocity and pressure values in the variables u_n and p_n respectively.
Now, if we want to calculate the stress values along the bottom horizontal boundary of the unit square. I have hardcoded the required values as follows
total_points = 10
x_coordinates = np.linspace(0,1,total_points)
y_coordinates = np.zeros(len(x_coordinates))
(Or else one can use, locate_dofs_topological or locate_dofs_geometrical):
Hence, IF I have the shear stress and normal stress vectors, then I can extract those values by:
domain = mesh
bb_tree = geometry.BoundingBoxTree(domain, domain.topology.dim)
points = np.zeros((3, len(x_coordinates)))
points[0] = x_coordinates
points[1] = y_coordinates
u_values = []
cells = []
points_on_proc = []
# Find cells whose bounding-box collide with the the points
cell_candidates = geometry.compute_collisions(bb_tree, points.T)
# Choose one of the cells that contains the point
colliding_cells = geometry.compute_colliding_cells(domain, cell_candidates, points.T)
for i, point in enumerate(points.T):
if len(colliding_cells.links(i))>0:
points_on_proc.append(point)
cells.append(colliding_cells.links(i)[0])
points_on_proc = np.array(points_on_proc, dtype=np.float64)
shear_stress_values = shear_stress.eval(points_on_proc, cells)
normal_stress_values = normal_stress.eval(points_on_proc, cells)
So, my task is to collect all the stress values to shear_stress and normal_stress
For this purpose, I construct function spaces as follows:
# Stresses with the variational form:
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())
n = FacetNormal(mesh)
T = -sigma(u_n, p_n)*n
Tn = inner(T, n) # Normal stress as a scalar
Tt = T - Tn*n # Shear stress as a vector
p_cg1 = FiniteElement('CG', mesh.ufl_cell(), 1)
scalar = FunctionSpace(mesh, p_cg1)
v_cg2 = VectorElement('CG', mesh.ufl_cell(), 2)
vector = FunctionSpace(mesh, v_cg2)
normal_stress = Function(scalar)
shear_stress = Function(vector)
Now to project the stress values to these function spaces, I was directed to this post: Obtain velocity normal to boundary
They have the following code written with FEniCS (not FEniCSx):
n = FacetNormal(mesh)
# Interpolate mesh values into boundary function
u_boundary = Function(D)
u_boundary.set_allow_extrapolation(True)
# Approximate facet normal in a suitable space using projection
n = FacetNormal(mesh)
V = VectorFunctionSpace(mesh, "CG", 2)
u_ = TrialFunction(V)
v_ = TestFunction(V)
a = inner(u_,v_)*ds
l = inner(n, v_)*ds
A = assemble(a, keep_diagonal=True)
L = assemble(l)
A.ident_zeros()
nh = Function(V)
solve(A, nh.vector(), L)
class normal_u(UserExpression):
def __init__(self, **kwargs):
super().__init__(**kwargs)
def eval(self, values, x):
n_eval = nh(x)
u_eval = u(x)
un = (u_eval[0]*n_eval[0] + u_eval[1]*n_eval[1])
values[0] = un * n_eval[0]
values[1] = un * n_eval[1]
def value_shape(self):
return (2,)
u_boundary.interpolate(normal_u())
So, by looking at it, I managed to write it upto:
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())
n = FacetNormal(mesh)
T = -sigma(u_n, p_n)*n
Tn = inner(T, n)
Tt = T - Tn*n
p_cg1 = FiniteElement('CG', mesh.ufl_cell(), 1)
scalar = FunctionSpace(mesh, p_cg1)
v_cg2 = VectorElement('CG', mesh.ufl_cell(), 2)
vector = FunctionSpace(mesh, v_cg2)
normal_stress = Function(scalar)
shear_stress = Function(vector)
u_ = TrialFunction(vector)
v_ = TestFunction(vector)
a = form(inner(u_,v_)*ds)
A = assemble_matrix(a)
A.assemble()
l = form(inner(n, v_)*ds)
L = create_vector(l)
Can someone please help me to write down the remaining matrix solver and then the class with user expression.
I’m extremely thankful for your help.