def METHOD_NAME(bathymetry, bnd_id, strength=8000.):
    """Smooths bathymetry near open boundaries"""
    fs = bathymetry.function_space()
    mesh = fs.mesh()
    solution = Function(fs, name='bathymetry')
    diffusivity = Function(fs, name='diff')
    # step 1: create diffusivity field
    delta_x = sqrt(CellVolume(mesh))
    distance = 2*delta_x
    get_boundary_relaxation_field(diffusivity, bnd_id, distance)
    # step 2: solve diffusion eq
    test = TestFunction(fs)
    f = inner(solution - bathymetry, test)*dx
    f += strength**2*diffusivity*inner(grad(solution), grad(test))*dx
    prob = NonlinearVariationalProblem(f, solution)
    solver = NonlinearVariationalSolver(prob)
    solver.solve()
    return solution