Parallel mesh refinement leads to corrupted solution when using read_checkpoint

General
1

Hi there.

I have the following problem. I’m trying to use write / read_checkpoint to load solutions back into Python in order to do some post-processing. This works fine in serial but in parallel my solutions return corrupted. After some debugging I found that my local mesh refinement function is causing the problem. The solutions can be loaded in Paraview fine but loading them back into Python using read_checkpoint leads to corruption. The refine() function seems to be the issue.

Any help?

Here is the main code:

from fenics import *

set_log_active(False)

a=1000

nx = 60
ny = 120
centrex = 0
centrey = 0.01*a
radius = 0.00229*a

x0 = 0
y0 = 0
x1 = 0.03*a
y1 = 0.06*a

def refiner(mesh, LEFT):
    cell_markers = MeshFunction("bool", mesh, mesh.topology().dim())
    cell_markers.set_all(False)
    for cell in cells(mesh):
        for facet in facets(cell): 
            for vertex in vertices(facet):
                if (vertex.point().array()[0] <= LEFT): 
                    cell_markers[cell] = True 
    mesh = refine(mesh, cell_markers)
    return mesh

def mesh_selector(mesh):
    mesh = refiner(mesh, 0.01*a)
    return mesh

mesh = RectangleMesh(Point(x0, y0), Point(x1, y1), nx, ny, 'left')
eps = Constant((2.5/6)*mesh.hmin())
mesh = mesh_selector(mesh)

dist = Expression('sqrt((pow((x[0]-A),2))+(pow((x[1]-B),2)))-r', degree=2, A=centrex, B=centrey, r=radius)
dist2 = Expression('(1/(1+exp((dist/eps))))', degree=2, eps=eps, dist=dist)

V = FunctionSpace(mesh, 'CG', 2)
phi0 = interpolate(dist2, V)

with XDMFFile("test/phi.xdmf") as outfile:
    outfile.write_checkpoint(phi0, "phi",0, append=False)

Here is the code I use to read in:

import matplotlib.pyplot as plt
fig = plt.figure(num=None, figsize=(8, 16), dpi=100,
                 facecolor='w', edgecolor='k')
ax1 = fig.add_subplot(1, 2, 1)
phi0=Function(V)
with XDMFFile('test/phi.xdmf') as infile:
    infile.read_checkpoint(phi0, "phi", 0)
plot(phi0)
ax1 = fig.add_subplot(1, 2, 2)
plot(interpolate(dist2, V))

Here is the result of that code:

output
output676×612 8.2 KB

2

I don’t have the bandwidth to try your code, but could your issue be related to this: Different result in serial and parallel run - #4 by nate ?

2 Likes
3

Hi Nate. This was indeed the solution I was looking for. Thank you very much.

1 Like