I’ve thought about this a bit more, and while the task I just posted is still a thing, my true problem is actually a bit more complicated…
Right now, I believe the code identifies the top 30% highest gradient magnitude cells in each rank, and then refines those cells. In this case, redistribution is not actually necessary, perhaps, as each rank will end up with the same kind of refinement.
However, the true problem that I, and I would think most folks doing this, would actually want to do, is to find the top 30% (or whatever fraction) cells across all ranks, and then refine those. In this case, there could be some MPI communication to view the entire global list of gradient magnitudes, get the top 30%, and then refine those cells within.
So, I guess now my task is to find out how I can use Allgather (or something) to analyze the global gradient magnitudes and their global cell indices, then pass that to relevant ranks, refine, and redistribute.
However, this does still lead me to the issue in the prior post, where refinement + redistribute leaves a mismatch between facet_indices and facet_markers.
EDIT:
Ok, I’ve solved this sub-task of refining based on a global list of gradient magnitudes, though it is probably not efficient (code below). There is probably a nice way to do it using the IndexMap functionality built into dolfinx, but I believe this works and gets back to the main task of getting meshtags for refined + redistributed meshes.
On this, I am still not sure what to do.
EDIT 2:
The parent_cell and parent_facet outputs seem to be the cells after partitioning, which means the ft values from the original gmshio.model_to_mesh are not helpful from this data.
It is looking like there’s a good reason I have seen no other code using gmsh meshtags with refinement + redistribute… is it even possible right now to do load-balancing mesh refinement while preserving meshtags?
#!/usr/bin/env python
from dolfinx import mesh
from dolfinx import fem
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import VTXWriter, gmshio
from dolfinx.fem import functionspace, dirichletbc, Function, Expression
from dolfinx.plot import vtk_mesh as create_vtk_mesh
import pyvista
import numpy as np
from petsc4py import PETSc
import ufl
from ufl import div, dx, grad, inner, dot
from mpi4py import MPI
from basix.ufl import element, mixed_element
import gmsh
comm = MPI.COMM_WORLD
if comm.rank == 0:
print('=========================================', flush=True)
print(' Setting up and solving new problem ', flush=True)
print('=========================================', flush=True)
# Size of mesh
H = 1.0
L = 1.0
# Domain dimension
gdim = 2
# Create mesh
gmsh.initialize()
gmsh.option.setNumber("General.Terminal",0)
if comm.rank == 0:
print('Generating mesh', flush=True)
# Create as rectangular domain with hole in the center
rectangle = gmsh.model.occ.addRectangle(0, 0, 0, L, H, tag=1)
# "Fluid" domain
fluid_marker = 1
gmsh.model.occ.synchronize()
volumes = gmsh.model.getEntities(dim=gdim)
assert (len(volumes) == 1)
gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker)
gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid")
# There will be an "inlet" on the left side, walls, and the cylindrical obstacle
inlet_marker, wall_marker = 2, 3
inflow, walls = [], []
boundaries = gmsh.model.getBoundary(volumes, oriented=False)
for boundary in boundaries:
center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
if np.allclose(center_of_mass, [0, H / 2, 0]):
inflow.append(boundary[1])
elif np.allclose(center_of_mass, [L / 2, H, 0]) or np.allclose(center_of_mass, [L / 2, 0, 0]) or np.allclose(center_of_mass, [L, H/2, 0]):
walls.append(boundary[1])
gmsh.model.addPhysicalGroup(1, walls, wall_marker)
gmsh.model.setPhysicalName(1, wall_marker, "Walls")
gmsh.model.addPhysicalGroup(1, inflow, inlet_marker)
gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
# Improve mesh near obstacle
res_max = H / 11
gmsh.option.setNumber("Mesh.MeshSizeMax", res_max)
gmsh.model.mesh.generate(gdim)
# converge gmsh model to doflinx mesh
msh, _, ft = gmshio.model_to_mesh(gmsh.model, comm, 0, gdim=gdim)
comm.Barrier()
# -------------------------------------------------------------------------
# ----------------------- COARSE FE SOLVE ---------------------------------
# -------------------------------------------------------------------------
# Set up and solve basic Poisson problem
if comm.rank == 0:
print('Setting up function space', flush=True)
V = functionspace(msh, ('Lagrange', 1))
# Dirichlet boundary condition at "inlet"
dofs = fem.locate_dofs_topological((V), msh.topology.dim-1, ft.find(2))
bc_inlet = fem.dirichletbc(PETSc.ScalarType(1), dofs, V)
# Walls
noslip = Function(V)
dofs = fem.locate_dofs_topological((V), msh.topology.dim-1, ft.find(3))
bc_walls = fem.dirichletbc(noslip, dofs)
bcs = [bc_inlet, bc_walls]
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = Function(V)
a = dot(grad(u), grad(v)) * dx
L = f * v * dx
# Solve
if comm.rank == 0:
print('Coarse solve', flush=True)
problem = LinearProblem(a, L, bcs = bcs, petsc_options={'ksp_type': 'preonly', 'pc_type': 'lu'})
uh = problem.solve()
# -------------------------------------------------------------------------
# --------------------- ADAPTIVE REFINEMENT -------------------------------
# -------------------------------------------------------------------------
if comm.rank == 0:
print('Adative refinement', flush=True)
# Mesh Refinement in FEniCSx
# Use current mesh level to locate cells worth refining
Vg = functionspace(msh, ("DG", 0))
F = Expression(ufl.sqrt(ufl.dot(ufl.grad(uh), ufl.grad(uh))),
Vg.element.interpolation_points()
)
grad_norm = Function(Vg)
grad_norm.interpolate(F)
# Local process gradient magnitudes
local_grad_mag = grad_norm.x.array
# BEGIN IDENTIFICATION OF GLOBAL TOP 30% GRADIENT MAG CELLS
# Task: identify the top 30% gradient magnitude cells across
# all ranks, and pass their local cell indices to a
# local mesh refinment
# Create a list filled with local process rank value
local_rank_list = np.full(local_grad_mag.shape, comm.rank)
# Create a list with local process cell indices
local_cell_list = list(range(local_grad_mag.shape[0]))
# Gather all gradient magnitudes, with their rank list
# and local cell index list
global_grad_mag = np.hstack(comm.allgather(local_grad_mag))
global_rank_list = np.hstack(comm.allgather(local_rank_list))
global_local_cell_list = np.hstack(comm.allgather(local_cell_list))
# Determine the top 30% gradient magnitude cell indices (global)
order = np.argsort(global_grad_mag)
global_cell_index = order[-int(0.3*order.size):-1]
# Determine which ranks these top 30% cells belong to
global_rank_list = global_rank_list[global_cell_index]
# Determine the local cell indices within those ranks
global_local_cell_index = global_local_cell_list[global_cell_index]
# For each local process, identify which (global) indices to
# pay attention to and refine
target_index = np.where(np.array(global_rank_list) == comm.rank)[0]
# Pick out the local cell index for those global indices
local_cell_index = global_local_cell_index[target_index]
# END GLOBAL TOP 30% IDENTIFICATION PROCESS
msh.topology.create_connectivity(1,2)
edges = mesh.compute_incident_entities(msh.topology, local_cell_index.astype(np.int32), 2, 1)
# gmsh refinement to keep facet markers
dim = msh.topology.dim
msh.topology.create_entities(dim-1)
msh.topology.create_entities(1)
meshtag = mesh.meshtags(
msh,
dim-1,
ft.indices,
ft.values,
)
if comm.rank == 0:
print('Refine plaza', flush=True)
fine_msh, parent_cell, parent_facet = mesh.refine_plaza(msh, edges, False, mesh.RefinementOption.parent_cell_and_facet)
fine_msh.topology.create_entities(dim-1)
if comm.rank == 0:
print('Transfer meshtags', flush=True)
new_ft = mesh.transfer_meshtag(meshtag, fine_msh, parent_cell, parent_facet)
if comm.rank == 0:
print('Create connectivity', flush=True)
fine_msh.topology.create_connectivity(fine_msh.topology.dim-1, fine_msh.topology.dim)
# Replacing old mesh with new refined mesh and meshtags for when this is
# done in a loop for multiple refinements
msh = fine_msh
ft = new_ft
# Visualize refined mesh
'''
if comm.rank == 0:
vtkdata = create_vtk_mesh(msh, msh.topology.dim)
plotter = pyvista.Plotter()
grid = pyvista.UnstructuredGrid(*vtkdata)
actor = plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
plotter.show()
'''
if comm.rank == 0:
print('Set up refined mesh function space', flush=True)
# Set up and solve basic Poisson problem
V = functionspace(msh, ('Lagrange', 1))
# Dirichlet boundary condition at "inlet"
dofs = fem.locate_dofs_topological((V), msh.topology.dim-1, ft.find(2))
bc_inlet = fem.dirichletbc(PETSc.ScalarType(1), dofs, V)
# Walls
noslip = Function(V)
dofs = fem.locate_dofs_topological((V), msh.topology.dim-1, ft.find(3))
bc_walls = fem.dirichletbc(noslip, dofs)
bcs = [bc_inlet, bc_walls]
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = Function(V)
a = dot(grad(u), grad(v)) * dx
L = f * v * dx
# Solve
if comm.rank == 0:
print('Solving refined problem', flush=True)
problem = LinearProblem(a, L, bcs = bcs, petsc_options={'ksp_type': 'preonly', 'pc_type': 'lu'})
uh = problem.solve()
# Write to file
data_fname = 'data2D_gmsh_refine_global.bp'
with VTXWriter(comm, data_fname, [uh]) as xf:
xf.write(0.0)
if comm.rank == 0:
print(f'Solution saved to file as {data_fname}', flush=True)
