Testfunctions integrate computes wrong results for Cylindrical2D coordinate system

[Da-Be-Ru]

It seems that the testfunction-based integration in VoronoiFVM supplies wrong results for the Cylindrical2D coordinate system. @chmerdon and I have provided a MWE with a hotfix that requires some discussion.

MWE

We calculate test functions $$T_i:\overline{\Omega}\rightarrow\mathbb{R}$$, $$i = 2,3$$ on $$\Omega=(0,2)\times (0,2)$$ such that $$T_i=1$$ on $$\Gamma_i$$ and $$T_i=0$$ on $$\partial\Omega\setminus \Gamma_i$$. We want to test the correct computation of the species mass transfer of a simple function $$u$$ across $$\Gamma_2=\{2\}\times [0,2]$$ or $$\Gamma_3=[0,2]\times\{2\}$$ both in the cartesian and the cylindrical case.

For a given source term $$s$$, we have that $$u$$ is a solution to the elliptic system

$$-\mathrm{div}(\mathbf{j}(u)(x,y)) = s(x,y) \text{ in } \, \Omega$$

in the cartesian case and

$$-\mathrm{div}_{\mathrm{cyl}}(\mathbf{j}_{\mathrm{cyl}}(u)(r,z)) = s(r,z) \text{ in } \, \Omega$$

in the axisymmetric cylindrical case (i.e. $$u$$ not varying in $$\theta$$) where $$\mathrm{div}_{\mathrm{cyl}}$$ the cylindrical divergence operator:

$$\mathrm{div}_{\mathrm{cyl}}(\mathbf{f})=\frac{1}{r}(\partial_r (r f_r) + \partial_z(r f_z))$$ $$\mathbf{j}_{\mathrm{cyl}}(u) = -\nabla_{(r,z)}(u) =- \partial_r u - \partial_z u.$$

With the flux $$\mathbf{j}(u)=-\nabla u$$ approximated via central differences, we expect for a simple function $$u_1(x,y) = x$$ or $$u_2(x,y)=y$$ (resp. $$u(r,z)=r$$ and $$u(r,z)=z$$ in the cylindrical case) the following:

Cartesian Case

$$\int_{\Gamma_i} \mathbf{j}(u_j)\cdot\mathbf{n}\mathrm{dS} = \begin{cases} \int_{\Gamma_2} -\nabla u_1 \cdot \mathbf{n} \, \mathrm{dS} = -2 & \text{ if } i=2, j=1 \\\int_{\Gamma_2} -\nabla u_2 \cdot \mathbf{n} \, \mathrm{dS} = 0 & \text{ if } i=2, j=2\\\int_{\Gamma_3}- \nabla u_1 \cdot \mathbf{n} \, \mathrm{dS} = 0 & \text{ if } i=3, j=1\\\int_{\Gamma_3} -\nabla u_2 \cdot \mathbf{n} \, \mathrm{dS} = -2 & \text{ if } i=3, j=2\end{cases}$$

Cylindrical Case

We denote by $$\hat{\Omega}=\Phi(\Omega\times [0,2\pi))$$ where $$\Phi(r,z,\theta) = [r\cos(\theta), r\sin(\theta),z]$$. In that case, we have

$$\begin{align} \int_{\hat{\Gamma_i}} \hat{\mathbf{j}}(u_j)\cdot\,\hat{\mathbf{n}}\,\mathrm{dS} &= 2\pi \int_{\Gamma_i} r\, \mathbf{j}_{\mathrm{cyl}}(u_j)\cdot \mathbf{n} \,\mathrm{dS} = \begin{cases}2\pi \int_{\Gamma_2} - r \nabla_{(r,z)} u_1 \cdot \mathbf{n} \, \mathrm{dS} = -8\pi & \text{ if } i=2, j=1 \\2\pi\int_{\Gamma_2} -r\nabla_{(r,z)} u_2 \cdot \mathbf{n} \, \mathrm{dS} = 0 & \text{ if } i=2, j=2\\2\pi\int_{\Gamma_3}- r\nabla_{(r,z)} u_1 \cdot \mathbf{n} \, \mathrm{dS} = 0 & \text{ if } i=3, j=1\\2\pi\int_{\Gamma_3} -r\nabla_{(r,z)} u_2 \cdot \mathbf{n} \, \mathrm{dS} = -4\pi & \text{ if } i=3, j=2\end{cases} \end{align}$$

Results

The Cartesian2D cases run through and deliver the expected results. However, in the Cylindrical2D case, the result for $$u_1$$ on $$\Gamma_2$$ and that of $$u_2$$ on $$\Gamma_3$$ are incorrect.

Possible Cause and Fix

It seems to us that there might be a missing corrective term for the mass flux in the cylindrical case. To circumvent this issue, we tried a hotfix by setting $$\tilde{T_i}(r,z) = rT_i(r,z)$$ and hijacking the integrate function for the Cartesian2D case:

$$\begin{align} \int_{\hat{\Gamma_i}} \hat{\mathbf{j}}(u_j)\cdot\,\hat{\mathbf{n}}\,\mathrm{dS} &= \int_{\hat{\Gamma}} \hat{T}_i\hat{\mathbf{j}}(u_j)\cdot\hat{\mathbf{n}}\,\mathrm{dS}\\\ &= 2\pi \int_{\Omega} \mathrm{div}_{\mathrm{cyl}}\left(T_i \mathbf{j}_{\mathrm{cyl}}(u_j) \right)\,r\,\mathrm{dr}\mathrm{dz} \\\ &= \int_{\Omega} \nabla_{(r,z)}\left(r\,T_i\right) \cdot \mathbf{j}_{\mathrm{cyl}}(u_j)\mathrm{dr}\mathrm{dz} + \int_{\Omega} rT_i \,\mathrm{div}_{(r,z)}(\mathbf{j}_{\mathrm{cyl}}(u_j)) \,\mathrm{dr}\mathrm{dz} \\\ &= \int_{\Omega} \nabla_{(r,z)}\tilde{T}_i \cdot \mathbf{j}_{\mathrm{cyl}}(u_j)\mathrm{dr}\mathrm{dz} + \int_{\Omega} \tilde{T}_i \,\mathrm{div}_{(r,z)}(\mathbf{j}_{\mathrm{cyl}}(u_j)) \,\mathrm{dr}\mathrm{dz} \end{align}.$$

Since the last term evaluates to zero, the first one in the sum would correspond to a call to integrate with a scaled test function $$\tilde{T} = rT$$.

This now produces a correct results. The only thing we keep wondering is that in the case of $$u_1$$ on $$\Gamma_3$$, we obtain a much worse accuracy than in the other cases where we expect $$0$$ as a result.

We included a test script with the hotfix, but we think this needs some discussion on how to fix this internally.

[Da-Be-Ru]

It seems we forgot to define a source term. VoronoiFVM assumes that $$u_i$$ is the solution to a corresponding elliptic problem whose terms are summed into the second summand in the last integral.

Will push an update on that one soon.

[Da-Be-Ru]

I now put in the correct source terms. The "corrective term" I spoke of earlier is not necessary - I simply made an error in the calculation.

The modified test is now set up thusly:

For each coordinate system, we have four different test cases, thus a total of 8 test cases.
For the Cartesian2D coordinate system, we have:

• For $$\Gamma_2$$, two functions $$c_{\Gamma_2,1}(x,y)=x$$ and $$c_{\Gamma_2,0}(x,y)=y$$ such that

$$\int_{\Gamma_2} \mathbf{j}(c_{\Gamma_2,1})\cdot\mathbf{n}\mathrm{dS} = \left|\Gamma_2\right| = 2 \quad \text{ and } \quad \int_{\Gamma_2}\mathbf{j}(c_{\Gamma_2,0})\cdot\mathbf{n}\mathrm{dS} = 0.$$

• For $$\Gamma_3$$, two functions $$c_{\Gamma_3,1}(x,y)=y$$ and $$c_{\Gamma_3,0}(x,y)=x$$ such that

$$\int_{\Gamma_3} \mathbf{j}(c_{\Gamma_3,1})\cdot\mathbf{n}\mathrm{dS} = \left|\Gamma_3\right| = 2 \quad \text{ and } \quad \int_{\Gamma_3}\mathbf{j}(c_{\Gamma_3,0})\cdot\mathbf{n}\mathrm{dS} = 0.$$

where $$\mathtt{source}(x,y)=-\mathrm{div}(\mathbf{j}(c_{\Gamma_{1,2},0,1})=0$$.
The integrate function approximates these as expected and everything is fine there.

For the Cylindrical2D coordinate system we have:

• For $$\Gamma_2$$, two functions $$c_{\Gamma_2,1}(r,z)=\frac{1}{2} \, r^2 - z^2$$ and $$c_{\Gamma_2,0}(r,z)=z$$ such that

$$2\pi \int_{\Gamma_2} r \mathbf{j}\_{\mathrm{cyl}}(c\_{\Gamma_2,1})\cdot\mathbf{n}\mathrm{dS} = 2\pi \int_0^2 2^2\,\mathrm{dz} = 16\pi \quad \text{ and } \quad \int_{\Gamma_2}\mathbf{j}(c_{\Gamma_2,0})\cdot\mathbf{n}\mathrm{dS} = 0,$$

where

$$\mathrm{div}_{\mathrm{cyl}}(\mathbf{j}_{\mathrm{cyl}}(c_{\Gamma_2,1}))(r,z)=1-2+\frac{1}{r} r = 0 \quad \text{ and } \quad \mathrm{div}_{\mathrm{cyl}}(\mathbf{j}_{\mathrm{cyl}}(c_{\Gamma_2,0}))(r,z)=0.$$

• For $$\Gamma_3$$, two functions $$c_{\Gamma_3,1}(r,z)=\frac{1}{2} \, r^2 - z^2$$ and $$c_{\Gamma_3,0}(r,z)=r$$ such that

$$2\pi \int_{\Gamma_3} r \mathbf{j}_{\mathrm{cyl}}(c_{\Gamma_3,1})\cdot\mathbf{n}\mathrm{dS} = 2\pi \int_0^2 r (-2\cdot2)\,\mathrm{dr} = -16\pi \quad \text{ and } \quad \int_{\Gamma_3}\mathbf{j}(c_{\Gamma_3,0})\cdot\mathbf{n}\mathrm{dS} = 0.$$

where

$$\mathrm{div}_{\mathrm{cyl}}(\mathbf{j}_{\mathrm{cyl}}(c_{\Gamma_3,1}))(r,z)=1-2+\frac{1}{r} r = 0 \quad \text{ and } \quad \mathrm{div}_{\mathrm{cyl}}(\mathbf{j}_{\mathrm{cyl}}(c_{\Gamma_3,0}))(r,z)=\frac{1}{r}.$$

Now the last example gives trouble since the source term contains $$-1/r$$ which turns into -Inf and is thus eventually passed into the mass flux from integrate.
I tried to address this by:

• disturbing $$r\leftarrow r+\varepsilon$$ for $$\varepsilon\approx 10^{-8}$$.
• replacing $$c_{\Gamma_3,0}(r,z)=\frac{rz^2}{2} - 2rz$$ which should also result in $$\int_{\Gamma_3}\mathbf{j}(c_{\Gamma_3,0})\cdot\mathbf{n}\mathrm{dS} = 0$$, but with a source term $$\mathrm{div}_{\mathrm{cyl}}(\mathbf{j}_{\mathrm{cyl}}(c_{\Gamma_3,0}))(r,z)=-r+\frac{z^2}{2r}-\frac{2z}{r}$$.

But neither of these yielded the desired result of $$0$$. I'm not quite sure if I've made a mistake here or if we even need such a comprehensive test. But I do think that these simple integrals should probably be evaluated correctly.