And here is a reproducible example with an built-in mesh.
You can vary the n_size and see that it works for, e.g. n_size=20 but there is a (small) error for n_size=200.
import dolfinx
from dolfinx.fem import FunctionSpace, Function, Constant, assemble_scalar, form, Expression
from dolfinx.io import XDMFFile
from dolfinx.mesh import locate_entities, meshtags, locate_entities_boundary
import dolfinx_mpc
import basix.ufl
import numpy as np
from mpi4py import MPI
from ufl import grad, dot, FacetNormal, Measure, TrialFunction, TestFunction, lhs, rhs, SpatialCoordinate
from petsc4py.PETSc import ScalarType
# Initialize Parallelization
############################
comm = MPI.COMM_WORLD
rank = comm.rank
def print0(string: str):
"""Print on rank 0 only"""
if MPI.COMM_WORLD.rank == 0:
print(string)
# mesh
######
n_size = 200
domain = dolfinx.mesh.create_unit_cube(
MPI.COMM_WORLD, n_size, n_size, n_size, cell_type=dolfinx.mesh.CellType.hexahedron)
radius = 0.1
midpoint = 0.5
def cells_sphere(x):
return ((x[0]-midpoint)**2 + (x[1]-midpoint)**2 + (x[2]-midpoint)**2) <= radius**2
def cells_cube(x):
return ((x[0]-midpoint)**2 + (x[1]-midpoint)**2 + (x[2]-midpoint)**2) >= radius**2
cells_0 = locate_entities(domain, domain.topology.dim, cells_sphere)
cells_1 = locate_entities(domain, domain.topology.dim, cells_cube)
# material parameters
#####################
S = FunctionSpace(domain, ("DG", 0))
kappa = Function(S)
kappa.x.array[:] = 2.5
kappa.x.array[cells_0] = np.full_like(cells_0, 5.0, dtype=ScalarType)
kappa.x.array[cells_1] = np.full_like(cells_1, 1.0, dtype=ScalarType)
# meshing measures
dx = Measure('dx', domain=domain) # Different volumes
x = SpatialCoordinate(domain)
domain.topology.create_connectivity(domain.topology.dim - 1,
domain.topology.dim)
# Element formulation
P1 = basix.ufl.element("Lagrange",
domain.topology.cell_name(), 1, gdim=domain.geometry.dim)
DG0 = basix.ufl.element("DG", domain.topology.cell_name(), 0, gdim=domain.geometry.dim, shape=(domain.geometry.dim,))
# Define FunctionSpaces
V = dolfinx.fem.functionspace(domain, P1)
W = dolfinx.fem.functionspace(domain, DG0)
T_ges = Function(V)
T_grad = Function(W)
# define trial and test functions
T_fluc = TrialFunction(V)
dT_fluc = TestFunction(V)
# boundary conditions
#####################
# prescribed "loads"
T_gradienten = []
for iii in ([1e-4]):
T_gradienten.append((iii, 0, 0))
T_gradienten.append((0, iii, 0))
T_gradienten.append((0, 0, iii))
T_gradienten = tuple(T_gradienten)
T_gradient = Constant(domain, ScalarType((0,0,0))) # Definition of a temperature gradient: can be changed later
bcs = []
# initialize Periodic Boundary Conditions
#########################################
pbc_directions = []
pbc_slave_tags = []
pbc_is_slave = []
pbc_is_master = []
pbc_meshtags = []
pbc_slave_to_master_maps = []
mpc = dolfinx_mpc.MultiPointConstraint(V)
def generate_pbc_slave_to_master_map(i):
def pbc_slave_to_master_map(x):
out_x = x.copy()
out_x[i] = x[i] - domain.geometry.x[:, i].max()
return out_x
return pbc_slave_to_master_map
def generate_pbc_is_slave(i):
return lambda x: np.isclose(x[i], domain.geometry.x[:, i].max())
def generate_pbc_is_master(i):
return lambda x: np.isclose(x[i], domain.geometry.x[:, i].min())
for i in range(domain.topology.dim):
pbc_directions.append(i)
pbc_slave_tags.append(i + 3)
pbc_is_slave.append(generate_pbc_is_slave(i))
pbc_is_master.append(generate_pbc_is_master(i))
pbc_slave_to_master_maps.append(generate_pbc_slave_to_master_map(i))
facets = locate_entities_boundary(domain, domain.topology.dim -1, pbc_is_slave[-1])
arg_sort = np.argsort(facets)
pbc_meshtags.append(meshtags(domain,
domain.topology.dim -1,
facets[arg_sort],
np.full(len(facets), pbc_slave_tags[-1], dtype=np.int32)))
N_pbc = len(pbc_directions)
for i in range(N_pbc):
# slave/master mapping of opposite surfaces (without slave-slave[-slave] intersections)
def pbc_slave_to_master_map(x):
out_x = pbc_slave_to_master_maps[i](x)
# remove edges that are connected to another slave surface
idx = np.logical_or(pbc_is_slave[(i + 1) % N_pbc](x),pbc_is_slave[(i + 2) % N_pbc](x))
out_x[i][idx] = np.nan
return out_x
mpc.create_periodic_constraint_topological(V, pbc_meshtags[i],
pbc_slave_tags[i],
pbc_slave_to_master_map,
bcs)
if len(pbc_directions) > 1:
def generate_pbc_slave_to_master_map_edges(dir_i, dir_j):
def pbc_slave_to_master_map_edges(x):
'''
Maps the slave edge dofs [intersection of slave_i, slave_j] (i=1,j=1) to
master corner dof [master_i, master_j] (i=0,j=0)
(1) map slave_x, slave_y to master_x, master_y: i=0, j=1
(2) map slave_x, slave_z to master_x, master_z: i=0, j=2
(3) map slave_y, slave_z to master_y, master_z: i=1, j=2
'''
out_x = x.copy()
out_x[dir_i] = x[dir_i] - domain.geometry.x[:, i].max()
out_x[dir_j] = x[dir_j] - domain.geometry.x[:, i].max()
idx = np.logical_and(pbc_is_slave[dir_i](x), pbc_is_slave[dir_j](x))
out_x[dir_i][~idx] = np.nan
out_x[dir_j][~idx] = np.nan
# remove corner DOFs at point (1,1,1)
idx_corner = np.logical_and(pbc_is_slave[0](x), np.logical_and(pbc_is_slave[1](x), pbc_is_slave[2](x)))
out_x[dir_i][idx_corner] = np.nan
out_x[dir_j][idx_corner] = np.nan
return out_x
return pbc_slave_to_master_map_edges
def pbc_slave_to_master_map_corner(x):
'''
Maps the slave corner dof [intersection of slave_x, slave_y, slave_z] (1,1,1) to
master corner dof [master_x, master_y, master_z] (0,0,0)
'''
out_x = x.copy()
out_x[0] = x[0] - domain.geometry.x[:, i].max()
out_x[1] = x[1] - domain.geometry.x[:, i].max()
out_x[2] = x[2] - domain.geometry.x[:, i].max()
idx_corner = np.logical_and(pbc_is_slave[0](x), np.logical_and(pbc_is_slave[1](x), pbc_is_slave[2](x)))
out_x[0][~idx_corner] = np.nan
out_x[1][~idx_corner] = np.nan
out_x[2][~idx_corner] = np.nan
return out_x
mapping_slave_intersections = [(0,1),(0,2),(1,2)] # pairs of slave intersections
# mapping slave intersection node (corner) to master intersection node (opposite corner)
mpc.create_periodic_constraint_topological(V, pbc_meshtags[0],
pbc_slave_tags[0],
pbc_slave_to_master_map_corner,
bcs)
for inters in mapping_slave_intersections:
# mapping slave intersections to opposite master intersections
mpc.create_periodic_constraint_topological(V, pbc_meshtags[inters[0]],
pbc_slave_tags[inters[0]],
generate_pbc_slave_to_master_map_edges(inters[0], inters[1]),
bcs)
mpc.finalize()
# weak formulation
##################
eqn = - kappa * dot(grad(dot(T_gradient , x) + T_fluc), grad(dT_fluc)) * dx
# solve
#######
a_form = lhs(eqn)
L_form = rhs(eqn)
petsc_options = {"ksp_rtol": 1e-9, "ksp_atol": 1e-9, "pc_type": "jacobi", "ksp_type": "cg"}
problem = dolfinx_mpc.LinearProblem(a_form, L_form, mpc, bcs=bcs, petsc_options=petsc_options)
i_n = 0
# get total volume
V_tot = domain.comm.allreduce(assemble_scalar(form(1 * dx)), op=MPI.SUM)
# create XDMF files
file_results_1 = XDMFFile(domain.comm, f"T_fluc_TestPerRB_{n_size}.xdmf", "w")
file_results_1.write_mesh(domain)
for T_grad_i in range(len(T_gradienten)):
i_n += 1
T_gradient.value = T_gradienten[T_grad_i]
T_fluc = problem.solve()
T_fluc.name = 'T_fluc'
# total temperature
T_ges_expr = Expression(dot(T_gradient, x) + T_fluc, V.element.interpolation_points())
T_ges.interpolate(T_ges_expr)
# temperature gradient
T_grad_expr = Expression(grad(T_ges), W.element.interpolation_points())
T_grad.interpolate(T_grad_expr)
# save to files
file_results_1.write_function(T_fluc, i_n)
# compute average temperature gradient
T_grad_avg = np.zeros(domain.topology.dim)
for k in range(domain.topology.dim):
T_grad_avg[k] = 1/V_tot * domain.comm.allreduce(assemble_scalar(form( T_grad[k] * dx)), op=MPI.SUM)
print0(f"Averaged temperature gradient in x: {T_grad_avg[0]} = {T_gradient.value[0]}?; in y: {T_grad_avg[1]} = {T_gradient.value[1]}?; in z: {T_grad_avg[2]} = {T_gradient.value[2]}? ")
file_results_1.close()
