Hi,everyone.
I am solving a nonlinear discontinuous finite element problem.
My speed field will gradually increase over time.I found during debugging that the solution of my first time step does not meet the boundary conditions.
In order for everyone to reproduce the problem, I am attaching my python code
import numpy as np
import ufl
from dolfinx import cpp as _cpp
from dolfinx import fem
from dolfinx.fem import (Constant, Function, FunctionSpace, dirichletbc,
extract_function_spaces, form,
locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.io import XDMFFile
from dolfinx.mesh import (CellType, GhostMode, create_rectangle,
locate_entities_boundary,create_unit_cube)
from ufl import div, dx, grad, inner,dot,outer,ds,dS,avg,jump,rot,cross
from mpi4py import MPI
from petsc4py import PETSc
# Create mesh
msh = create_unit_cube(MPI.COMM_WORLD,
8,8,8,
CellType.tetrahedron, GhostMode.none)
u_element=ufl.FiniteElement("RT", ufl.tetrahedron, 1) # H（div） Vh
A_element=ufl.FiniteElement("N1curl",ufl.tetrahedron,1) # H(curl) Ch
P_element=ufl.FiniteElement("DG",ufl.tetrahedron,0) # Qh
TH= ufl.MixedElement([u_element, A_element,P_element])
W =fem.FunctionSpace(msh,TH)
W0,_=W.sub(0).collapse()
W1,_=W.sub(1).collapse()
W2,_=W.sub(2).collapse()
(un,An,Pn)=ufl.TrialFunctions(W) #
(v,fai,q)=ufl.TestFunctions(W) #
#Vh=fem.FunctionSpace(msh,u_element)
un_1=fem.Function(W0)
un_2=fem.Function(W0)
def initial_u(x):
values=np.zeros((3,x.shape[1]))
values[0,:]=x[1] #y
values[1,:]=x[2] #z
values[2,:]=x[0] #x
return values
un_1.interpolate(initial_u)
#Ch=fem.FunctionSpace(msh,A_element)
An_1=fem.Function(W1)
An_2=fem.Function(W1)
def initial_A(x):
values=np.zeros((3,x.shape[1]))
values[0,:]=x[2] #z
values[1,:]=x[0]*0 #0
values[2,:]=x[1] #y
return values
An_1.interpolate(initial_A)
P=fem.Function(W2)
def initial_P(x):
values=np.zeros((x.shape[1]))
return values
P.interpolate(initial_P)
#---------------------------------------------------------------------------------------
x = ufl.SpatialCoordinate(msh)
h = ufl.CellDiameter(msh)
n = ufl.FacetNormal(msh)
unv=(un+un_1)*0.5
uns=1.5*un_1-0.5*un_2
Anv=(An+An_1)*0.5
Ans=1.5*An_1-0.5*An_2
lmbda = ufl.conditional(ufl.gt(ufl.dot(uns, n), 0), 1, 0)
u_uw = lmbda("+") * unv("+") + lmbda("-") * unv("-")
tao_1=1000
k=1
Rm=1
alph=10
b = 2*tao_1*inner(unv,v)*dx
b -= inner(unv,div(outer(uns,v)))*dx
b += inner((dot(uns, n))("+") * u_uw, v("+")) * dS+inner((dot(uns, n))("-") * u_uw, v("-")) * dS #
b += inner(dot(uns, n) * lmbda * unv, v) * ds
b += inner(grad(unv),grad(v))*dx +alph/h('+')*inner(jump(unv),jump(v))*dS
b -= inner(avg(dot(grad(unv),n)),jump(v))*dS-inner(avg(dot(grad(v),n)),jump(unv))*dS
b += 2*k*tao_1*inner(rot(Anv),rot(fai))*dx
b += k*inner(2*tao_1*Anv+cross(rot(Ans),unv),2*tao_1*fai+cross(rot(Ans),v))*dx
b += inner(2*tao_1*un_1,v)*dx+k*inner(2*tao_1*An_1,2*tao_1*fai+cross(rot(Ans),v))*dx
b += 2*tao_1*inner(div(unv),div(v))*dx+inner(div(unv),q)*dx-inner(div(v),Pn)*dx
a=ufl.lhs(b)
L=ufl.rhs(b)
a=fem.form(a)
L=fem.form(L)
#---------------------------------------------------------------------------------
def local_A(x):
return np.isclose(x[1],0) #mark y=1
def local_u(x):
return np.logical_and(np.isclose(x[2],0),np.isclose(x[1],0)) #mark z=0 and y=0
def local_p(x):
return np.logical_and(np.isclose(x[2],0),np.isclose(x[1],0),np.isclose(x[0],0)) #mark z=0 and x=0
def u_D(x):
values=np.zeros((3,x.shape[1]))
values[0,:]=0 #y
values[1,:]=0 #z
values[2,:]=x[0] #x
return values
def A_D(x):
values=np.zeros((3,x.shape[1]))
values[0,:]=x[2] #z
values[1,:]=0 #z
values[2,:]=0 #x
return values
def p_D(x):
values=np.zeros((1,x.shape[1]))
values[0,:]=0 #z
return values
Vh=fem.FunctionSpace(msh,u_element)
u_boundary = Function(W0)
u_boundary.interpolate(u_D)
facets = locate_entities_boundary(msh,2, local_u)
dofs = locate_dofs_topological((W.sub(0), W0), 2, facets)
bc0 = dirichletbc(u_boundary, dofs, W.sub(0))
Dh=fem.FunctionSpace(msh,A_element)
A_boundary = Function(W1)
A_boundary.interpolate(A_D)
facets = locate_entities_boundary(msh,2, local_A)
dofs = locate_dofs_topological((W.sub(1), W1), 2, facets)
bc1 = dirichletbc(A_boundary, dofs, W.sub(1))
Ph=fem.FunctionSpace(msh,P_element)
p_boundary = Function(W2)
p_boundary.interpolate(p_D)
facets = locate_entities_boundary(msh,2, local_u)
dofs = locate_dofs_topological((W.sub(2), W2), 2, facets)
bc2 = dirichletbc(p_boundary, dofs, W.sub(2))
bcs=[bc0,bc1,bc2]
#--------------------------------------------------------------------------------------
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
create_vector, set_bc)
A=assemble_matrix(a,bcs=bcs)
A.assemble()
b=create_vector(L)
'''solver = PETSc.KSP().create(msh.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.BCGS)
pc1 = solver.getPC()
print(pc1)
pc1.setType(PETSc.PC.Type.HYPRE)
pc1.setHYPREType("boomeramg")
solver.setFromOptions()'''
solver = PETSc.KSP().create(msh.comm)
solver.setOperators(A)
solver.setType("preonly")
solver.getPC().setType("lu")
solver.getPC().setFactorSolverType("superlu_dist")
opts = PETSc.Options() # type: ignore
opts["ksp_error_if_not_converged"] = 1
def domain_average(msh, v):
"""Compute the average of a function over the domain"""
vol = msh.comm.allreduce(fem.assemble_scalar(fem.form(fem.Constant(msh, 1.0) * dx)), op=MPI.SUM)
return (1 / vol) * msh.comm.allreduce(fem.assemble_scalar(fem.form(v * dx)), op=MPI.SUM)
wn=Function(W)
t=0
for n in range(1):
if n==0:
t += 0.001
un_2.x.array[::]=un_1.x.array[::]
An_2.x.array[::]=An_1.x.array[::]
fem.petsc.assemble_matrix(a,bcs=bcs) # type: ignore
assemble_vector(b, L)
apply_lifting(b, [a], [bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, bcs)
solver.solve(b, wn.vector)
u = wn.sub(0).collapse()
A = wn.sub(1).collapse()
p = wn.sub(2).collapse()
un_2.x.array[::]=un_1.x.array[::]
un_1.x.array[::]=u.x.array[::]
An_2.x.array[::]=An_1.x.array[::]
An_1.x.array[::]=A.x.array[::]
un_2.x.scatter_forward()
un_1.x.scatter_forward()
An_2.x.scatter_forward()
An_1.x.scatter_forward()
p.x.array[:] -= domain_average(msh, p)
else:
t += 0.001
fem.petsc.assemble_matrix(a,bcs=bcs) # type: ignore
with b.localForm() as b_loc:
b_loc.set(0)
assemble_vector(b, L)
apply_lifting(b, [a], [bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, bcs)
solver.solve(b, wn.vector)
u = wn.sub(0).collapse()
A = wn.sub(1).collapse()
p = wn.sub(2).collapse()
un_2.x.array[::]=un_1.x.array[::]
un_1.x.array[::]=u.x.array[::]
An_2.x.array[::]=An_1.x.array[::]
An_1.x.array[::]=A.x.array[::]
un_2.x.scatter_forward()
un_1.x.scatter_forward()
An_2.x.scatter_forward()
An_1.x.scatter_forward()
p.x.array[:] -= domain_average(msh, p)
print('t is',t)
u = wn.sub(0).collapse()
A = wn.sub(1).collapse()
p = wn.sub(2).collapse()
with XDMFFile(MPI.COMM_WORLD, "out_WH/velocity.xdmf", "w") as ufile_xdmf:
u.x.scatter_forward()
ufile_xdmf.write_mesh(msh)
ufile_xdmf.write_function(u)
with XDMFFile(MPI.COMM_WORLD, "out_WH/A.xdmf", "w") as pfile_xdmf:
A.x.scatter_forward()
pfile_xdmf.write_mesh(msh)
pfile_xdmf.write_function(A)
with XDMFFile(MPI.COMM_WORLD, "out_WH/p.xdmf", "w") as pfile_xdmf:
p.x.scatter_forward()
pfile_xdmf.write_mesh(msh)
pfile_xdmf.write_function(p)
The following are the node coordinates that match the boundary：
The rows are x, y, and z respectively.
The following are the corresponding velocity field values
The rows are u_x, u_y, and u_z respectively.
This question has been bothering me for a long time, and any answer will be appreciated 