Wrong results in ver-1a case?

[RemDelaporteMathurin]

Hi all!

I recently looked at the ver-1a verification case and I noticed the solution differs from that of the TMAP7 V&V.

This is the fractional release of TMAP8:

This is the fractional release of TMAP7:

I've looked at the parameters and they all look the same in TMAP8 and TMAP7 (page C-3 of V&V document):

• diffusivity
• solubility
• temperature
• geometry (enclosure volume, surface and length of slab)
• initial pressure

Could you comment on this difference? Is there a parameter that I missed?

Thanks in advance!!

[singhgp4321]

@RemDelaporteMathurin The case is taken from TMAP4. Is TMAP7 ver-1a same as TMAP4 ver-1a?

[RemDelaporteMathurin]

Sorry I may have gotten mixed up in links. I'll recap:

TMAP4:

TMAP7:

TMAP8

It is clear here that the TMAP8 results are from TMAP4 and not TMAP7 as I mentionned. Apologies.

Note: i think the $T$ in the exponential here should be lowercase $t$.

Now, when I try to compare the TMAP8 results to the analytical solution (equation) given in the docs and in TMAP4 ref, it doesn't produce the plotted results. See below:

See code

```python import numpy as np import matplotlib.pyplot as plt k = 1.38065e-23 # J/mol Boltzmann constant def get_roots_bis(L, alpha_max, step=0.0001): """Gets the roots of alpha = L / tan(alpha) Args: L (float): parameter L alpha_max (float): the maximum alpha to consider step (float, optional): the step discretising alphas. The smaller the step, the more accurate the roots. Defaults to 0.0001. Returns: np.array: array of roots """ alphas = np.arange(0, alpha_max, step=step)[1:] f = alphas g = L / np.tan(alphas) plt.plot(alphas, f, "-") plt.plot(alphas, g, "-") idx = np.argwhere(np.diff(np.sign(f - g))).flatten() # remove one every other idx idx = idx[::2] plt.plot(alphas[idx], f[idx], "ro") plt.show() roots = alphas[idx] return roots # FIXME this doesn't work # taken from https://mooseframework.inl.gov/TMAP8/verification/ver-1a.html#longhurst1992verification def analytical_expression_FR(t, D, S, V, T, A, l): """ Args: t (float, ndarray): time (s) D (float): diffusivity (m2/s) S (float): solubility (H/m3/Pa) V (float): enclosure volume (m3) T (float): temperature (K) A (float): enclosure surface area (m2) l (float): slab length (m) """ phi = 1 / (k * T * S) L = l * A / (V * phi) roots = get_roots_bis(L=L, alpha_max=2000, step=3e-4) roots = roots[:, np.newaxis] sec = 1 / np.cos(roots) summation = (2 * L * sec - np.exp(-(roots**2) * D * t / l**2)) / ( L * (L + 1) + roots**2 ) summation = np.sum(summation, axis=0) FR = 1 - summation return FR if __name__ == "__main__": T = 2373 times = np.linspace(0, 45, 1000) FR_tmap8 = analytical_expression_FR( t=times, D=2.6237e-11, S=7.244e22 / T, V=5.20e-11, T=T, A=2.16e-6, l=3.3e-5, ) plt.plot(times, FR_tmap8, label="TMAP8") plt.ylim(bottom=0) plt.ylabel("Fractional release") plt.xlabel("Time (s)") plt.show() ```

However, when plotting the cumulative release from the enclosure (ie the integral of the particle flux at x=0) I obtain the correct curve:

See code

```python import numpy as np from scipy.integrate import cumtrapz import matplotlib.pyplot as plt # https://stackoverflow.com/questions/28766692/how-to-find-the-intersection-of-two-graphs/28766902#28766902 k = 1.38065e-23 # J/mol Boltzmann constant def get_roots(L, l, alpha_max, step=0.0001): """Gets the roots of alpha = L / tan(alpha * l) Args: L (float): parameter L l (float): parameter l alpha_max (float): the maximum alpha to consider step (float, optional): the step discretising alphas. The smaller the step, the more accurate the roots. Defaults to 0.0001. Returns: np.array: array of roots """ alphas = np.arange(0, alpha_max, step=step)[1:] f = alphas g = L / np.tan(alphas * l) # plt.plot(alphas, f, "-") # plt.plot(alphas, g, "-") idx = np.argwhere(np.diff(np.sign(f - g))).flatten() # remove one every other idx idx = idx[::2] # plt.plot(alphas[idx], f[idx], "ro") # plt.show() roots = alphas[idx] return roots def analytical_expression_flux(t, P_0, D, S, V, T, A, l): """ value of the flux at the external surface (not in contact with pressure) J = -D * dc/dx Args: t (float, ndarray): time (s) P_0 (float): initial presure (Pa) D (float): diffusivity (m2/s) S (float): solubility (H/m3/Pa) V (float): enclosure volume (m3) T (float): temperature (K) A (float): enclosure surface area (m2) l (float): slab length (m) """ L = S * T * A * k / V roots = get_roots(L=L, l=l, alpha_max=1e7, step=1) roots = roots[:, np.newaxis] summation = (np.exp(-(roots**2) * D * t) * roots) / ( (l * (roots**2 + L**2) + L) * np.sin(roots * l) ) last_term = summation[-1] summation = np.sum(summation, axis=0) print(last_term / summation) flux = 2 * S * P_0 * L * D * summation return flux def cumulative_flux(t, P_0, D, S, V, T, A, l): flux = analytical_expression_flux(t, P_0, D, S, V, T, A, l) cumulative_flux = cumtrapz(flux, t, initial=0) initial_quantity = P_0 * V / k / T return cumulative_flux * A / initial_quantity if __name__ == "__main__": T = 2373 times = np.linspace(0, 45, 1000) cum_flux = cumulative_flux( t=times, P_0=1e6, D=2.6237e-11, S=7.244e22 / T, V=5.20e-11, T=T, A=2.16e-6, l=3.3e-5, ) plt.plot(times, cum_flux) plt.ylabel("cumulative flux") plt.xlabel("Time (s)") plt.show() ```

My question is therefore:

Could you please provide the script with which you obtain the analytical solution values in the analytical.csv file (assuming it's not simply copy pasted from the TMAP4 report)?

Would it be possible that the authors of the TMAP7 report notice an error in the analytical expression of the TMAP4 report which would have justified changing the metric that is used to verify the code?

Thank you so much for your support on this - very - intriguing issue.

[RemDelaporteMathurin]

Update:

I tried to derive the expression for the cumulative flux shown above.

Starting from the analytical solution of the concentration given in the TMAP7 report (Equation 1):

$$ c(x,t) = 2 L P{0} S \sum{n=1}^{\infty} \frac{e^{- D t \alpha_n^2} \sin{\left(x \alpha_n \right)}}{\left(L + l \left(L^{2} + \alpha_n^2\right)\right) \sin{\left(l \alpha_n \right)}} $$

where $\alpha_n$ are the roots of

$$ \alpha_n = \frac{L}{\tan(\alpha_n \ l)} $$

NOTE: the $\alpha_n$ are different in this report compared to TMAP4/TMAP8

The cumulative release $\mathrm{R}$ is given by:

$$ \mathrm{R} = \int_0^t J(x=0, t) dt = \int0^t D \nabla c |{x=0} dt $$

From the expression of the concentration we can obtain:

$$ \mathrm{R}= 2 L P{0} S \sum{n=1}^{\infty} \frac{1 - e^{- D t_{f} \alpha_n^{2}}}{\left(L + l \left(L^{2} + \alpha_n^2\right)\right) \alpha_n\sin{\left(l \alpha_n \right)}} $$

Now multiplying by the surface area and normalising by the initial quantity of particles in the enclosure, the fractional release $\mathrm{FR}$ is:

$$ \mathrm{FR} = \frac{\mathrm{R} \ A}{P_0 \ V / (k_B T)} $$

Finally:

$$ \mathrm{FR} = \frac{2 L S A}{V / (kB T)} \sum{n=1}^{\infty} \frac{1 - e^{- D t_{f} \alpha_n^{2}}}{\left(L + l \left(L^{2} + \alpha_n^2\right)\right) \alpha_n\sin{\left(l \alpha_n \right)}} $$

I plotted this expression against the one of the TMAP4/TMAP8 docs as well as the expression from TMAP7.

If I'm not mistaken this expression is very different from the one given in the TMAP4/TMAP8. And yet it seems to reproduce the analytical results obtained by TMAP4/TMAP8. The code to produce this is below.

See code

```python import numpy as np import matplotlib.pyplot as plt # https://stackoverflow.com/questions/28766692/how-to-find-the-intersection-of-two-graphs/28766902#28766902 k = 1.38065e-23 # J/mol Boltzmann constant def get_roots(L, l, alpha_max, step=0.0001): """Gets the roots of alpha = L / tan(alpha * l) Args: L (float): parameter L l (float): parameter l alpha_max (float): the maximum alpha to consider step (float, optional): the step discretising alphas. The smaller the step, the more accurate the roots. Defaults to 0.0001. Returns: np.array: array of roots """ alphas = np.arange(0, alpha_max, step=step)[1:] f = alphas g = L / np.tan(alphas * l) # plt.plot(alphas, f, "-") # plt.plot(alphas, g, "-") idx = np.argwhere(np.diff(np.sign(f - g))).flatten() # remove one every other idx idx = idx[::2] # plt.plot(alphas[idx], f[idx], "ro") # plt.show() roots = alphas[idx] return roots def get_roots_bis(L, alpha_max, step=0.0001): """Gets the roots of alpha = L / tan(alpha) Args: L (float): parameter L alpha_max (float): the maximum alpha to consider step (float, optional): the step discretising alphas. The smaller the step, the more accurate the roots. Defaults to 0.0001. Returns: np.array: array of roots """ alphas = np.arange(0, alpha_max, step=step)[1:] f = alphas g = L / np.tan(alphas) plt.plot(alphas, f, "-") plt.plot(alphas, g, "-") idx = np.argwhere(np.diff(np.sign(f - g))).flatten() # remove one every other idx idx = idx[::2] plt.plot(alphas[idx], f[idx], "ro") plt.show() roots = alphas[idx] return roots def analytical_expression_fractional_release_TMAP7(t, P_0, D, S, V, T, A, l): """ FR = 1 - P(t) / P_0 where P(t) is the pressure at time t and P_0 is the initial pressure Taken from https://inldigitallibrary.inl.gov/sites/sti/sti/4215153.pdf Equation 4 Note: in the report, the expression of FR is given as P(T)/P_0, but it shown as 1 - P(t)/P_0 in the graph (Figure 1) Args: t (float, ndarray): time (s) P_0 (float): initial presure (Pa) D (float): diffusivity (m2/s) S (float): solubility (H/m3/Pa) V (float): enclosure volume (m3) T (float): temperature (K) A (float): enclosure surface area (m2) l (float): slab length (m) """ L = S * T * A * k / V roots = get_roots(L=L, l=l, alpha_max=1e6, step=1) print(len(roots)) roots = roots[:, np.newaxis] summation = np.exp(-(roots**2) * D * t) / (l * (roots**2 + L**2) + L) last_term = summation[-1] summation = np.sum(summation, axis=0) print(last_term / summation) pressure = 2 * P_0 * L * summation fractional_release = 1 - pressure / P_0 return fractional_release # FIXME this doesn't work # taken from https://mooseframework.inl.gov/TMAP8/verification/ver-1a.html#longhurst1992verification def analytical_expression_fractional_release_TMAP8(t, D, S, V, T, A, l): """ Analytical expression for the fractional release given by TMAP4 report and TMAP8 docs This doesn't produce the correct analytical results presented by TMAP4 and TMAP8. Note, this expression isn't used in the TMAP7 V&V report. see https://github.com/idaholab/TMAP8/issues/75 for discussion with TMAP8 team Args: t (float, ndarray): time (s) D (float): diffusivity (m2/s) S (float): solubility (H/m3/Pa) V (float): enclosure volume (m3) T (float): temperature (K) A (float): enclosure surface area (m2) l (float): slab length (m) """ phi = 1 / (k * T * S) L = l * A / (V * phi) roots = get_roots_bis(L=L, alpha_max=2000, step=3e-4) roots = roots[:, np.newaxis] sec = 1 / np.cos(roots) summation = (2 * L * sec - np.exp(-(roots**2) * D * t / l**2)) / ( L * (L + 1) + roots**2 ) last_term = summation[-1] summation = np.sum(summation, axis=0) print(summation[0]) print(last_term / summation) fractional_release = 1 - summation return fractional_release def analytical_expression_flux(t, P_0, D, S, V, T, A, l): """ value of the flux at the external surface (not in contact with pressure) J = -D * dc/dx Args: t (float, ndarray): time (s) P_0 (float): initial presure (Pa) D (float): diffusivity (m2/s) S (float): solubility (H/m3/Pa) V (float): enclosure volume (m3) T (float): temperature (K) A (float): enclosure surface area (m2) l (float): slab length (m) """ L = S * T * A * k / V roots = get_roots(L=L, l=l, alpha_max=1e7, step=1) roots = roots[:, np.newaxis] summation = (np.exp(-(roots**2) * D * t) * roots) / ( (l * (roots**2 + L**2) + L) * np.sin(roots * l) ) last_term = summation[-1] summation = np.sum(summation, axis=0) print(last_term / summation) flux = 2 * S * P_0 * L * D * summation return flux def cumulative_flux_analytical(t, P_0, D, S, V, T, A, l): """ Analytical expression for the cumulative flux at the external surface (not in contact with pressure) integral(A D grad(c) . n dt) between 0 and t normalised by the initial quantity in the enclosure Args: t (float, ndarray): time (s) P_0 (float): initial presure (Pa) D (float): diffusivity (m2/s) S (float): solubility (H/m3/Pa) V (float): enclosure volume (m3) T (float): temperature (K) A (float): enclosure surface area (m2) l (float): slab length (m) """ L = S * T * A * k / V roots = get_roots(L=L, l=l, alpha_max=1e7, step=1) roots = roots[:, np.newaxis] num = 1 - np.exp(-(roots**2) * D * t) denum = roots * np.sin(roots * l) * (l * (roots**2 + L**2) + L) summation = num / denum last_term = summation[-1] summation = np.sum(summation, axis=0) print(last_term / summation) cumulative_flux_val = 2 * S * P_0 * L * summation initial_quantity = P_0 * V / k / T return cumulative_flux_val * A / initial_quantity if __name__ == "__main__": T = 2373 times = np.linspace(0, 45, 1000) cum_flux = cumulative_flux_analytical( t=times, P_0=1e6, D=2.6237e-11, S=7.244e22 / T, V=5.20e-11, T=T, A=2.16e-6, l=3.3e-5, ) FR_tmap8 = analytical_expression_fractional_release_TMAP8( t=times, D=2.6237e-11, S=7.244e22 / T, V=5.20e-11, T=T, A=2.16e-6, l=3.3e-5, ) FR_tmap7 = analytical_expression_fractional_release_TMAP7( t=times, P_0=1e6, D=2.6237e-11, S=7.244e22 / T, V=5.20e-11, T=T, A=2.16e-6, l=3.3e-5, ) FR_TMAP8 = np.genfromtxt("analytical.csv", delimiter=",", names=True) plt.plot(times, cum_flux, label="new formula", linewidth=3, color="tab:green") plt.plot(times, FR_tmap7, label="TMAP7") plt.plot(times, FR_tmap8, label="TMAP4/TMAP8 (equation)", color="tab:red") plt.scatter( FR_TMAP8["times"], FR_TMAP8["frac_rel"], label="TMAP8 (csv)", color="tab:green", alpha=0.5, ) plt.ylabel("FR") plt.xlabel("Time (s)") plt.legend() plt.show() ```

Maybe an error was made in the derivation of the TMAP4 analytical solution? It's impossible to tell as the authors of the report give no details (or references) on how they obtained this analytical solution.

[RemDelaporteMathurin]

Hi! Any update on this? @cticenhour

[cticenhour]

Hi Remi! I am unaware of whether @singhgp4321 has had the chance to look at your explorations on these discrepancies yet, but I have set aside time in the morning to take a look at this and will report back when I have something useful to contribute. 😄

[singhgp4321]

@RemDelaporteMathurin I plotted the analytical curve for TMAP4 results based on the numbers they have in the report, and it is not based on their analytical expression. I am not sure why the values from the analytical expression don't match their results . . . May be there is an error in the expression, as you indicate.

[RemDelaporteMathurin]

@singhgp4321 thanks for the reply

Whenever possible, especially in a V&V, one should avoid simply copy pasting things from past reports and hard-coding data like this. The result is that the documentation page is misleading to users.

Do you confirm that the figure of merit you are computing and comparing to the analytical solution is the integral of the cumulative flux on the non-exposed surface normalised by the initial quantity of H in the enclosure?

I believe it is very important to clarify this and fix any wrong equation to avoid confusing users like me! V&V should be a reliable and robust process. If one claims an certain equation is solved, then one should solve/implement this equation, not take results from decades ago for granted.

✅ Some verification cases do that correctly like the ver-1b case, the ver-1d case, the ver-1fa case, the ver-1fb case, and the ver-1g case

❌ But a couple of cases (representing almost 40% of the verification cases!) have a hard-coded analytical solution that could very well suffer the same issues like ver-1c: https://github.com/idaholab/TMAP8/blob/22acf97325e7fa040e6ac44741a52f4e5f88ebb2/test/tests/ver-1c/comparison_ver-1c.py#L20

ver-1e: https://github.com/idaholab/TMAP8/blob/22acf97325e7fa040e6ac44741a52f4e5f88ebb2/test/tests/ver-1e/comparison_ver-1e.py#L18-L20

[singhgp4321]

@RemDelaporteMathurin Are you sure that the analytical solution in TMAP4 is wrong?

[RemDelaporteMathurin]

Are you sure that the analytical solution in TMAP4 is wrong?

I mean you are more than welcome to try and plot it yourself. I shared all the code to reproduce the figures I showed. I'd love to be proven wrong - maybe I made a mistake somewhere.

That's why I asked if you guys had plotted the analytical solution given in TMAP4 - I tried, it didn't work. I tried rederiving the analytical solution myself - and it works. Nothing more I can do really.

[simopier]

@RemDelaporteMathurin, thank you for pointing these things out and being thorough in your description of the potential issue! I'll look into your comments and the cause of the differences between TMAP4/7 later this week and come back to you with an explanation or a fix.

[cticenhour]

@RemDelaporteMathurin To give a brief update on my end - yesterday, I did get to speak to @humrickhouse about the discrepancies here, and he suggested that the green and blue functions in your last plot (from TMAP4 and TMAP7 manuals respectively) may be representations of fractional release from different sides of the enclosure. He said he wasn't sure why the lead author for the TMAP7 manual decided to change the plotted quantity between releases, but that he may chime in here if he has time. I had less time than I thought yesterday, so will be looking through this today. As @simopier said, thanks again for pointing this out, and we'll get back to you soon!

[RemDelaporteMathurin]

the green and blue functions in your last plot (from TMAP4 and TMAP7 manuals respectively) may be representations of fractional release from different sides of the enclosure

That is possible. In the TMAP7 report, $x=0$ is the side not exposed to the gas where $c=0$ and $x=l$ is the surface exposed to the gas.

I'm more interested in comparing the red and green curves (and points). The way I rederived the analytical equation is computing the fractional release at $x=0$ (from the TMAP7 analytical solution of the concentration field). Looking at the TMAP8 file, the metric is the flux integral on the right hand side:

https://github.com/idaholab/TMAP8/blob/c711bc80a10bddee5319d0d560a10a63dd5346e3/test/tests/ver-1a/ver-1a.i#L84-L92

And the right hand side corresponds to $c=0$:

https://github.com/idaholab/TMAP8/blob/c711bc80a10bddee5319d0d560a10a63dd5346e3/test/tests/ver-1a/ver-1a.i#L42-L47

Which means that the equation given in TMAP4/TMAP8 should be equivalent to the equation I derived above. Again, the most straightforward way of unravelling the mystery here is to simply plot the equation given in TMAP4/TMAP8... 😄

[cticenhour]

Interesting. Agreed! 😄 Will chime back in soon.

[simopier]

I started to dig into ver-1a and created a PR with some changes and some discussions around the next steps to address this issue: #86.

[RemDelaporteMathurin]

As discussed in #86 there was a typo in my code. Both the analytical solution I derived from the TMAP7 solution and the TMAP4 solutions were equivalent!

Closing this.

[simopier]

I will add to what @RemDelaporteMathurin said that significant improvements have been done to the ver-1a case to include verification case 1a from TMAP7 (in addition to the existing 1a from TMAP4), increase transparency/documentation/reproducibility, and fix other typos as a response to this issue being created.

The PR https://github.com/idaholab/TMAP8/pull/86 with these changes is currently under review and will be pushed once it satisfies reviewers.

Thank you all (especially @RemDelaporteMathurin) again for your engagement and contribution.