Hi everyone !
When using discontinuous finite element method（DG） for interior surface integral(dS),we must use (‘-’) and (‘+’) to restrict.
So I want to know the principle of dolfinx partitioning (‘-’)and(‘+’)?
Further,When running programs in parallel,will (‘-’) and (‘+’) change after grid partitioning?
Thanks dokken!
Running an interior surface integral in parallel yields different results than running it in series.
Is this because after partitioning the grid, (‘-’) and (‘+’) have changed?
Are you using dolfinx.mesh.GhostMode.shared_facet when reading in/creating your mesh?
Usually in DG methods, there are a combination of avg and jump terms which means that the orientation doesn’t matter.
The consistent “+” and “-” orientation is usually only important if you try to do a one-sided interior integral (only one restriction), or if different things happens on different sides of an interface
I use dolfinx.mesh GhostMode.shared_facet
read_mesh = io.XDMFFile(MPI.COMM_WORLD,"kk.xdmf" ,"r")
msh=read_mesh.read_mesh(name="Grid")
It be used as the default parameter.
I need to use one-sided interior integral (dS) instead of ds.
The serial program, I have tested and it is correct.
If it is caused by grid partitioning, can I customize grid partitioning to solve it?
I would need a minimal reproducible example to see what goes wrong in your case, as it there shouldn’t be any parallel problems with this approach.
My grid is a spherical mesh, and the three spherical surfaces with radii ri(=0.538),ro(=1.538) and 5*ro divide the mesh into three subdomains.
This example only focuses on the second subdomain (ri<r<ro).
Firstly, I create custom interior and boundary surfaces tag：
ds should be used on the blue surface.
(1/Pr)* inner(grad(mean(T,T_n_1)),outer(Tv,n))* ds
Replace with
(1/Pr)* inner(grad(mean(T,T_n_1)),outer(Tv,n))(‘-’) * dS(2)
(1/Pr)* inner(grad(mean(T,T_n_1)),outer(Tv,n))(‘+’) * dS(3)
Serial operation is not a problem, and the result is also correct.
Parallel running iterative solver error.
But when I commented it out ：
F+=-Ekman * inner(grad(mean(u,u_n_1)), outer(v, n))('-') * dS(2)
F+=-Ekman * inner(grad(mean(u,u_n_1)), outer(v, n))('+') * dS(3)
F+=-Ekman * inner(outer(mean(u,u_n_1), n), grad(v))('-') * dS(2)
F+=-Ekman * inner(outer(mean(u,u_n_1), n), grad(v))('+') * dS(3)
F+=Ekman * (alpha / h('-')) * inner(outer(mean(u,u_n_1), n), outer(v, n))('-') * dS(2)
F+=Ekman * (alpha / h('+')) * inner(outer(mean(u,u_n_1), n), outer(v, n))('+') * dS(3)
The solver works normally.
doflinx version :0.7
use grid is:
https://easylink.cc/o7kwsc
https://easylink.cc/a9m131
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
from dolfinx import default_real_type, fem, io, mesh,default_scalar_type
from dolfinx.fem.petsc import assemble_matrix, assemble_vector
from ufl import (CellDiameter, FacetNormal, TestFunction, TrialFunction, avg,
conditional, div, dot, dS, ds, dx, grad, gt, inner, outer,
TrialFunctions,TestFunctions,curl,cross)
import ufl
from dolfinx.io import gmshio,XDMFFile
from basix.ufl import element,mixed_element
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
create_vector, set_bc,create_matrix,)
import time
import os
import dolfinx
def norm_div(comm, v):
return np.sqrt(comm.allreduce(fem.assemble_scalar(fem.form(inner(div(v), div(v)) * dx)), op=MPI.SUM))
def creat_measure(integral_type,msh,dim,entitiesmarkers):
entities_indices, entities_markers = [], []
for (marker, locator) in entitiesmarkers:
entities = mesh.locate_entities(msh, dim, locator)
entities_indices.append(entities)
entities_markers.append(np.full_like(entities, marker))
entities_indices = np.hstack(entities_indices).astype(np.int32)
entities_markers = np.hstack(entities_markers).astype(np.int32)
sorted_entities = np.argsort(entities_indices)
entities_tag = mesh.meshtags(msh, dim, entities_indices[sorted_entities], entities_markers[sorted_entities])
if dim==msh.topology.dim-1:
msh.topology.create_connectivity(msh.topology.dim-1, msh.topology.dim)
return ufl.Measure(integral_type, domain=msh, subdomain_data=entities_tag),entities_tag
def creat_tag(msh,dim,entitiesmarkers):
entities_indices, entities_markers = [], []
for (marker, locator) in entitiesmarkers:
entities = mesh.locate_entities(msh, dim, locator)
entities_indices.append(entities)
entities_markers.append(np.full_like(entities, marker))
entities_indices = np.hstack(entities_indices).astype(np.int32)
entities_markers = np.hstack(entities_markers).astype(np.int32)
sorted_entities = np.argsort(entities_indices)
entities_tag = mesh.meshtags(msh, dim, entities_indices[sorted_entities], entities_markers[sorted_entities])
if dim==msh.topology.dim-1:
msh.topology.create_connectivity(msh.topology.dim-1, msh.topology.dim)
return entities_tag
def norm_r(x):
return (x[0]**2+x[1]**2+x[2]**2)**0.5
def r_2(x):
return (x[0]**2+x[1]**2+x[2]**2)
def noslip(x):
return np.stack((np.zeros(x.shape[1]), np.zeros(x.shape[1]), np.zeros(x.shape[1])))
def boundT_0(x):
values = np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType)
return values
def boundT_1(x):
values = np.ones((1, x.shape[1]), dtype=PETSc.ScalarType)
return values
def T_x(x):
return (2*norm_r(x)-ri-ro)
def initial_T(x):
return ro*ri/norm_r(x)-ri+21/np.sqrt(17920*np.pi)*(1-3*T_x(x)**2+3*T_x(x)**4-T_x(x)**6)*((x[0]**2+x[1]**2)/(x[0]**2+x[1]**2+x[2]**2))**2*\
(1-8*x[0]**2/(x[0]**2+x[1]**2)+8*x[0]**4/(x[0]**2+x[1]**2)**2)
###########
num_time_steps =100
t_end = 0.1
Ekman=0.001
Ra=100
Pr=1
ro=20/13
ri=7/13
outside_radius_ratio=5
out_interval=1
read_mesh = io.XDMFFile(MPI.COMM_WORLD,"cp6n_0.45k_full.xdmf" ,"r")
msh=read_mesh.read_mesh(name="Grid")
gdim = msh.geometry.dim
fdim = msh.topology.dim - 1
domains=[(11,lambda x: x[0]**2+x[1]**2+x[2]**2>=ro**2-0.0001),
(12,lambda x: (x[0]**2+x[1]**2+x[2]**2<=ro**2+0.0001)&(x[0]**2+x[1]**2+x[2]**2>=ri**2-0.0001)),
(13,lambda x: x[0]**2+x[1]**2+x[2]**2<=ri**2+0.0001)]
boundaries = [(1, lambda x: np.isclose(x[0]**2+x[1]**2+x[2]**2, (outside_radius_ratio*ro)**2)),
(2, lambda x: np.isclose(x[0]**2+x[1]**2+x[2]**2, ro**2)),
(3, lambda x: np.isclose(x[0]**2+x[1]**2+x[2]**2, ri**2))]
interiors = [ (2, lambda x: np.isclose(x[0]**2+x[1]**2+x[2]**2, ro**2)),
(3, lambda x: np.isclose(x[0]**2+x[1]**2+x[2]**2, ri**2)),
(24,lambda x: (x[0]**2+x[1]**2+x[2]**2<=ro**2+0.001)&(x[0]**2+x[1]**2+x[2]**2>=ri**2-0.001))]
dx,cell_tag=creat_measure("dx",msh,gdim,domains)
ds,facet_tag=creat_measure("ds",msh,fdim,boundaries)
interiors_tag=creat_tag(msh,fdim,interiors)
dx_oc=dx(12)
indices1=interiors_tag.find(24)
indices2=interiors_tag.find(2)
indices3=interiors_tag.find(3)
indices23 = np.union1d(indices2, indices3)
indices_in = np.setdiff1d(indices1, indices23)
indices=[]
indices.append(indices_in)
indices.append(indices2)
indices.append(indices3)
mark=[]
mark.append(np.full_like(indices_in,33))
mark.append(np.full_like(indices2,2))
mark.append(np.full_like(indices3,3))
marks=np.hstack(mark).astype(np.int32)
indices_in=np.hstack(indices).astype(np.int32)
sorted_indices = np.argsort(indices_in)
msh.topology.create_connectivity(msh.topology.dim-1, msh.topology.dim)
tags_in_facets = mesh.meshtags(msh, fdim , indices_in[sorted_indices], marks[sorted_indices])
dS=ufl.Measure('dS', domain=msh, subdomain_data=tags_in_facets)
dS_in=dS(33)
with io.XDMFFile(msh.comm, "interiors_nb2.xdmf", "w") as xdmf:
xdmf.write_mesh(msh)
xdmf.write_meshtags(tags_in_facets, msh.geometry)
print('tag_in is finished')
U_element=element("RT", msh.basix_cell(), 2)
P_element=element("DG", msh.basix_cell(), 1)
T_element=element("DG", msh.basix_cell(), 1)
MIXED= mixed_element([U_element, P_element,T_element])
W =fem.functionspace(msh,MIXED)
W0,Wmap0=W.sub(0).collapse()
W1,Wmap1=W.sub(1).collapse()
W2,Wmap2=W.sub(2).collapse()
(u,P,T)=ufl.TrialFunctions(W)
(v,q,Tv)=ufl.TestFunctions(W)
imap = msh.topology.index_map(gdim)
num_cells = imap.size_local
ghost_cells = imap.num_ghosts
print(f"ID {msh.comm.rank}, num owned cells {num_cells} num ghost cells {ghost_cells},num dofs {W.dofmap.index_map.size_local * W.dofmap.index_map_bs}", flush=True)
class BoundaryCondition_u():
def __init__(self, type, marker, values):
self._type = type
if type == "Dirichlet":
u_D = fem.Function(W0)
u_D.interpolate(values)
facets = facet_tag.find(marker)
dofs = fem.locate_dofs_topological((W.sub(0),W0), fdim, facets)
self._bc = fem.dirichletbc(u_D, dofs,W.sub(0))
elif type == "Neumann":
self._bc = inner(values, v) * ds(marker)
elif type == "Robin":
self._bc = values[0] * inner(u-values[1], v)* ds(marker)
else:
raise TypeError("Unknown boundary condition: {0:s}".format(type))
@property
def bc(self):
return self._bc
@property
def type(self):
return self._type
class BoundaryCondition_T():
def __init__(self, type, marker, values):
self._type = type
if type == "Dirichlet":
u_D = fem.Function(W2)
u_D.interpolate(values)
dofs = fem.locate_dofs_geometrical((W.sub(2), W2), boundaries[marker-1][1])
# fem.locate_dofs_topological((W.sub(2),W2), fdim, facets)
self._bc = fem.dirichletbc(u_D, dofs,W.sub(2))
elif type == "Neumann":
self._bc = inner(values, v) * ds(marker)
elif type == "Robin":
self._bc = values[0] * inner(u-values[1], v)* ds(marker)
else:
raise TypeError("Unknown boundary condition: {0:s}".format(type))
@property
def bc(self):
return self._bc
@property
def type(self):
return self._type
boundary_conditions_u = [BoundaryCondition_u("Dirichlet", 2, noslip),
BoundaryCondition_u("Dirichlet", 3, noslip)]
boundary_conditions_T = [BoundaryCondition_T("Dirichlet", 2, boundT_0),
BoundaryCondition_T("Dirichlet", 3, boundT_1)]
bcs = []
for condition in boundary_conditions_u:
if condition.type == "Dirichlet":
bcs.append(condition.bc)
for condition in boundary_conditions_T:
if condition.type == "Dirichlet":
bcs.append(condition.bc)
u_n_1=fem.Function(W0)
u_n_2=fem.Function(W0)
u_D = fem.Function(W0)
A_n_1=fem.Function(W2)
A_n_2=fem.Function(W2)
T_n_1=fem.Function(W2)
T_n_2=fem.Function(W2)
T_n_1.interpolate(initial_T,cell_tag.find(12))
T_n_2.interpolate(initial_T,cell_tag.find(12))
T_D=fem.Function(W2)
T_D.interpolate(boundT_1)
T_D0=fem.Function(W2)
solution=fem.Function(W)
delta_t = fem.Constant(msh, default_real_type(t_end / num_time_steps))
alpha = fem.Constant(msh, default_real_type(6.0))
h = CellDiameter(msh)
n = FacetNormal(msh)
x = ufl.SpatialCoordinate(msh)
def mean(vn,vn_1):
return 1/2*(vn+vn_1)
def star(v_n_1,v_n_2):
return 1/2*(3*v_n_1-v_n_2)
def jump(phi, n):
return outer(phi("+"), n("+")) + outer(phi("-"), n("-"))
def u_devied_curl(u):
return ufl.as_vector([u[0]/(u[2].dx(1)-u[1].dx(2)),
u[1]/(u[0].dx(2)-u[2].dx(0)),
u[2]/(u[1].dx(0)-u[0].dx(1))])
lmbda = conditional(gt(dot(star(u_n_1,u_n_2), n), 0), 1, 0)
u_uw = lmbda("+") * mean(u,u_n_1)("+") + lmbda("-") * mean(u,u_n_1)("-")
T_uw = lmbda("+") * mean(T,T_n_1)("+") + lmbda("-") * mean(T,T_n_1)("-")
zeros=fem.Constant(msh, default_real_type(0))
### N-S equation
F=Ekman*(inner(u / delta_t, v) * dx_oc - inner(u_n_1 / delta_t, v) * dx_oc)
##trilinear form
F+= Ekman*-inner(mean(u,u_n_1), div(outer(v, star(u_n_1,u_n_2)))) * dx_oc
F+= Ekman*inner((dot(star(u_n_1,u_n_2), n))("+") * u_uw, v("+")) * dS_in
F+= Ekman*inner((dot(star(u_n_1,u_n_2), n))("-") * u_uw, v("-")) * dS_in
F+= Ekman*inner(dot(star(u_n_1,u_n_2), n) * lmbda * u, v)('-') * dS(2)
F+= Ekman*inner(dot(star(u_n_1,u_n_2), n) * lmbda * u, v)('+') * dS(3)
####Viscous
F+=Ekman * inner(grad(mean(u,u_n_1)), grad(v)) * dx_oc
F+=-Ekman * inner(avg(grad(mean(u,u_n_1))), jump(v, n)) * dS_in
F+=-Ekman * inner(jump(mean(u,u_n_1), n), avg(grad(v))) * dS_in
F+=Ekman * (alpha / avg(h)) * inner(jump(mean(u,u_n_1), n), jump(v, n)) * dS_in
F+=-Ekman * inner(grad(mean(u,u_n_1)), outer(v, n))('-') * dS(2)
F+=-Ekman * inner(grad(mean(u,u_n_1)), outer(v, n))('+') * dS(3)
F+=-Ekman * inner(outer(mean(u,u_n_1), n), grad(v))('-') * dS(2)
F+=-Ekman * inner(outer(mean(u,u_n_1), n), grad(v))('+') * dS(3)
F+=Ekman * (alpha / h('-')) * inner(outer(mean(u,u_n_1), n), outer(v, n))('-') * dS(2)
F+=Ekman * (alpha / h('+')) * inner(outer(mean(u,u_n_1), n), outer(v, n))('+') * dS(3)
F+=2*inner(cross(ufl.as_vector([0,0,1]),mean(u,u_n_1)),v)*dx_oc
F+=- inner(P, div(v)) * dx_oc - inner(div(u), q) * dx_oc
F -= Ra * inner(ufl.as_vector([x[0],x[1],x[2]])*star(T_n_1,T_n_2),v) *dx_oc
F += inner(T / delta_t, Tv) * dx_oc -inner(T_n_1 / delta_t, Tv) * dx_oc
F += -inner(mean(T,T_n_1), div(outer(Tv, star(u_n_1,u_n_2)))) * dx_oc
F+=inner((dot(star(u_n_1,u_n_2), n))("+") * T_uw, Tv("+")) * dS_in
F+=inner((dot(star(u_n_1,u_n_2), n))("-") * T_uw, Tv("-")) * dS_in
F+=inner(dot(star(u_n_1,u_n_2), n) * lmbda * T, Tv)('-') * dS(2)
F+=inner(dot(star(u_n_1,u_n_2), n) * lmbda * T, Tv)('+') * dS(3)
F+= (1/Pr)*inner(grad(mean(T,T_n_1)), grad(Tv)) * dx_oc
F+=- (1/Pr)* inner(avg(grad(mean(T,T_n_1))), jump(Tv, n)) * dS_in
F+=- (1/Pr)* inner(jump(mean(T,T_n_1), n), avg(grad(Tv))) * dS_in
F+= (1/Pr)* (alpha / avg(h)) * inner(jump(mean(T,T_n_1), n), jump(Tv, n)) * dS_in
F+=- (1/Pr)* inner(grad(mean(T,T_n_1)),outer(Tv,n))('-') * dS(2)
F+=- (1/Pr)* inner(grad(mean(T,T_n_1)),outer(Tv,n))('+') * dS(3)
F+=-(1/Pr)* inner(outer(mean(T,T_n_1), n), grad(Tv))('-') * dS(2)
F+=-(1/Pr)* inner(outer(mean(T,T_n_1), n), grad(Tv))('+') * dS(3)
F+= (1/Pr)* (alpha / h('-')) * inner(mean(T,T_n_1), Tv)('-') * dS(2)
F+= (1/Pr)* (alpha / h('+')) * inner(mean(T,T_n_1), Tv)('+') * dS(3)
F+=inner(zeros*u,v)*dx
F+=inner(zeros*P,q)*dx
F+=inner(zeros*T,Tv)*dx
aa=ufl.lhs(F)
LL=ufl.rhs(F)
a=fem.form(aa)
L=fem.form(LL)
A = create_matrix(a)
b=create_vector(L)
def make_reasons(reason):
return dict([(getattr(reason,r),r)
for r in dir(reason) if not r.startswith('_')])
def get_ksp_reason(ksp):
pc_re=make_reasons(PETSc.PC.FailedReason())
r_pc=ksp.getPC().getFailedReason()
return pc_re[r_pc]
def compute_ksp(ksp):
values=ksp.computeExtremeSingularValues()
return values
ksp = PETSc.KSP().create(msh.comm)
ksp.setOperators(A)
ksp.setType("gmres")
ksp.getPC().setType(PETSc.PC.Type.BJACOBI)
opts = PETSc.Options()
opts["ksp_rtol"] = 1.0e-10
opts["ksp_atol"] = 1.0e-3
opts["ksp_max_it"] =2000
ksp.setFromOptions()
def monitor_residual(ksp, its, rnorm):
if MPI.COMM_WORLD.Get_rank() ==0:
print(f'Iteration {its}, Residual norm {rnorm}', flush=True)
ksp.setMonitor(monitor_residual)
'''ksp = PETSc.KSP().create(msh.comm) # type: ignore
ksp.setOperators(A)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("mumps")
opts = PETSc.Options() # type: ignore
opts["mat_mumps_icntl_14"] = 80 # Increase MUMPS working memory
opts["mat_mumps_icntl_24"] = 1 # Option to support solving a singular matrix (pressure nullspace)
opts["mat_mumps_icntl_25"] = 0 # Option to support solving a singular matrix (pressure nullspace)
opts["ksp_error_if_not_converged"] = 1
ksp.setFromOptions()'''
A.zeroEntries()
fem.petsc.assemble_matrix(A, a, bcs=bcs)
A.assemble()
zerosrow=A.findZeroRows()
x = A.createVecRight()
A.zeroRows(zerosrow.array[::],diag=1.0,b=x)
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)
zerosrow=A.findZeroRows()
print('zero rows is',len(zerosrow.array[::]))
ksp.solve(b, solution.vector)
u_h=solution.sub(0).collapse()
p_h=solution.sub(1).collapse()
T_h=solution.sub(2).collapse()
e_u=norm_div(msh.comm,u_h)
print('div(u) is',e_u)
print('pc_error type',get_ksp_reason(ksp))