Critique on "Volume penalization for inhomogeneous Neumann boundary conditions modeling scalar flux in complicated geometry"
In this letter, we provide counter-examples to demonstrate that it is possible to retain second-order accuracy using Sakurai et al.'s method, even when different flux boundary conditions are imposed on multiple interfaces that do not conform to the Cartesian grid. We consider both continuous and discontinuous indicator functions in our test problems. Both indicator functions yield a similar convergence rate for the problems considered here. We also find that the order of accuracy results for some of the cases presented in Sakurai et al. are not reproducible. This is demonstrated by re-considering the same one- and two-dimensional Poisson problems solved in Sakurai et al. in this letter. The results shown in this letter demonstrate that the spatial order of accuracy of the flux-based VP approach of Sakurai et al. is between $\mathcal{O}$(1) and $\mathcal{O}$(2), and it depends on the underlying problem/model. The spatial order of accuracy cannot simply be deduced a priori based on the imposed flux values, shapes, or grid-conformity of the interfaces, as concluded in Sakurai et al. Further analysis is required to understand the spatial convergence rate of the flux-based VP method.