Hello everyone.I am trying to solve an MHD problem：（Actually, it’s example 6.1 from
A fully divergence-free finite element method formagnetohydrodynamic equations
DOI: 10.1142/S0218202518500173）
The variational form is:
where
the python code is :
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
#fai_element=ufl.FiniteElement("DG",ufl.tetrahedron,2) # SH
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(un_1, n), 0), 1, 0)
u_uw = lmbda("+") * un_1("+") + lmbda("-") * un_1("-") # u -> upwind flux
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(un_1, n))("+") * u_uw, v("+")) * dS+inner((dot(un_1, n))("-") * u_uw, v("-")) * dS #
b += inner(dot(un_1, n) * lmbda * un, 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
a=ufl.lhs(b)
L=ufl.rhs(b)
a=fem.form(a)
L=fem.form(L)
#------------------------------------------------------------------------------------
def local_A(x): #mark boundary of A
return np.isclose(x[1],0) #mark y=1
def local_u(x): #mark boundary of u
return np.logical_and(np.isclose(x[2],0),np.isclose(x[0],0)) #mark z=0 && 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
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), Vh), 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), Dh), 2, facets)
bc1 = dirichletbc(A_boundary, dofs, W.sub(1))
bcs=[bc0,bc1]
#----------------------------------------------------------------------------------------------
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()
pc1.setType(PETSc.PC.Type.HYPRE)
pc1.setHYPREType("boomeramg")
solver.setFromOptions()
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)
I want to know why does the solution I calculate not even meet the Dirichlet boundary condition?
The analytical solution is：
