Add ver-1kc-2

[lkadz]

(Ref. #12)

[moosebuild]

Job Documentation, step Sync to remote on aea9aea wanted to post the following:

View the site here

This comment will be updated on new commits.

[simopier]

This is not a full review, I'm just providing suggestions to solve the current test failures.

Also, make sure to add a reference to the issue number in the PR message at the top of the page.

[lkadz]

Hi @simopier ,

I had a question regarding the reaction $H_2 + T_2 \longleftrightarrow 2HT$. In TMAP8, using MatReactionFlexible, we rely on a kinetic law (e.g. for HT, $\frac{d(c_{HT})}{dt} = 2K \cdot c_{H_2} \cdot c_{T_2}$), whereas in TMAP7, we use an equilibrium formula: $P_{HT} = \eta \cdot \bigl(P_{H_2}\bigr)^{0.5} \cdot \bigl(P_{T_2}\bigr)^{0.5}$

Could you please advise on how to ensure the equilibrium law is respected in TMAP8? I tried to illustrate this in the last Python subplot figure, but the equilibrium isn’t fully reached yet.

Also, I set an arbitrary reaction rate constant K and would appreciate any insight into how its value is usually determined. I compared it to other cases that use MatReactionFlexible but couldn’t figure out how those values were obtained.

Thanks in advance for your help!
Best

[moosebuild]

Job Coverage, step Generate coverage on aea9aea wanted to post the following:

Coverage

Coverage did not change

Full coverage report

This comment will be updated on new commits.

[simopier]

You're on the right track!
Here's a partial review with some suggestions. I'll need to do a deeper review once you have all the pieces together.

[simopier]

Hi @simopier ,

I had a question regarding the reaction H 2 + T 2 ⟷ 2 H T . In TMAP8, using MatReactionFlexible, we rely on a kinetic law (e.g. for HT, d ( c H T ) d t = 2 K ⋅ c H 2 ⋅ c T 2 ), whereas in TMAP7, we use an equilibrium formula: P H T = η ⋅ ( P H 2 ) 0.5 ⋅ ( P T 2 ) 0.5

Could you please advise on how to ensure the equilibrium law is respected in TMAP8? I tried to illustrate this in the last Python subplot figure, but the equilibrium isn’t fully reached yet.

Also, I set an arbitrary reaction rate constant K and would appreciate any insight into how its value is usually determined. I compared it to other cases that use MatReactionFlexible but couldn’t figure out how those values were obtained.

Great questions!

You correctly identified that while TMAP7 uses an equilibrium condition, we probably want to use a kinetic description in TMAP8 to be more general. However, we need to reconcile the two to obtain the same equilibrium condition.
You currently struggle to reconcile the kinetic equation with the equilibrium formula because you are missing a piece in your kinetic equation. The reaction goes in both directions, and you are currently missing the decomposition of HT into H2 and T2:
$\frac{dc_{HT}}{dt} = 2 K_1 c_{H2}c_{T2} - K_2c_{HT}^2$
and similarly:
$\frac{dc_{H2}}{dt} = - K_1 c_{H2}c_{T2} + \frac{1}{2} K_2c_{HT}^2$
and
$\frac{dc_{T2}}{dt} = - K_1 c_{H2}c_{T2} + \frac{1}{2} K_2c_{HT}^2$

at equilibrium:
$\frac{dc_{HT}}{dt} = \frac{dc_{HT}}{dt} = \frac{dc_{HT}}{dt} = 0$,
which leads to
$2 K_1 c_{H2}c_{T2} - K_2c_{HT}^2 = 0$,
which can be rearranged into
$P_{HT}=\eta \sqrt{P_{H2}P_{T2}}$
with
$\eta = \sqrt{\frac{2K_1}{K_2}}$.

So, to obtain the same equilibrium as in TMAP7, you can define $K_1$ and $K_2$ in several ways, as long as you get $\eta = \sqrt{\frac{2K_1}{K_2}}$. However, if you also want to ensure that the equilibrium condition is obtained quickly knowing that it is assumed to be immediate in TMAP7, then you'll want to define $K_1$ and $K_2$ large enough so that it is much faster than other kinetics such as diffusion and/or surface reactions or other kinetics relevant to whatever case you are working on at that time.

Please double check the math above, but that should give you all you need to address your questions above. Here, you'll need to add blocks for the chemical reactions (feel free to check other cases if you are unsure how to do it), and select appropriate values for $K_1$ and $K_2$.

I hope that helps! Let me know if you have other questions.

[lkadz]

Hi @simopier,
Thank you for your valuable insights! In the updated commits, I have ensured that the chemical reactions respect the TMAP7 equilibrium.
In this branch (ver-1kc-2), I renamed ver-1kc.md to ver-1kc-1.md, but I noticed there is still a documentation error. Should I also add ver-1kc-1.md to the ver-1kc branch to resolve this issue?

[simopier]

The ver-1kc branch has already been merged in your previous PR. Right now, you are getting failures because you have not updated ver-1kc.md into ver-1kc-1.md everywhere. For example, in the test specification file for ver-1kc-1, it still says

verification = `ver-1kc.md`

instead of

verification = `ver-1kc-1.md`

[simopier]

Here's a new round of reviews. Let me know if you have any questions.

[simopier]

The mass variation shown in the figure below is actually much larger, and the results do not correspond to the TMAP7 predictions. The kernels you have defined might not conserve mass, probably because one uses the wrong coefficients or added/missed a term - or maybe your formula for mas conservation is incorrect.
Make sure to resolve this issue, explain any difference with the TMAP7 case if some still remain, and update this line of the documentation with the latest number.

[simopier]

Things are looking good!

My main remaining concerns are:

• the very loose tolerances. They should be much smaller.
• the fact that we do not match the right equilibrium condition (very likely due to large tolerances and large time steps)
• The discussion about K1 and K2 in the documentation, where we are not exactly on the same page.

As you try to refine the tolerances, feel free to play with the values of K_1 and K_2. Try to make them lower and higher (still with eta=2) to understand their impact on the results and on convergence.

[simopier]

Suggested change

	dt = 0.01
	dt = 1e-3

[simopier]

Suggested change

	growth_factor = 1.05
	growth_factor = 1.1

[simopier]

Suggested change

	cutback_factor_at_failure = 0.8
	cutback_factor = 0.9
	cutback_factor_at_failure = 0.9

[simopier]

These tolerances are very large. Try to get them closer to what we commonly use.

[simopier]

Suggested change

	dtmax = 10
	dtmax = 0.1

[simopier]

Suggested change

	simulation_time = '${units 1 s}'
	simulation_time = '${units 0.25 s}'

Not much seems to happen after that.

[simopier]

The equilibrium is much better with tighter tolerances! But I think we still have some room for improvement, especially for enclosure 2.

What was your experience playing with the tolerances and seeing their impact on the results?

[simopier]

Suggested change

	should_execute = false # this test relies on the output files from ver-1kc-2_csv, so it shouldn't be run twice
	should_execute = false # this test relies on the output files from ver-1kc-2_csv_heavy, so it shouldn't be run twice

[simopier]

Why are you using such large numbers for these two?
This means that for fewer than 20-10=10 iterations, the timestep is increased, nothing is done between 20-10=10 and 20+10=30 iterations, and the timestep is refined above 20+10=30 iterations.
However, note that since you have nl_max_its = 20, it actually never goes above 20, and immediately refines the timestep above 20 iterations.

These numbers are usually set up lower than this, especially the iteration window.

[lkadz]

A significant number of nonlinear convergences occur for iterations exceeding 10. By setting nl_max_its = 20 and preventing the decrease in time steps between 10 and 20 (with optimal_iterations = 20 and iteration_window = 10), I observed that the time steps increased sufficiently to achieve quicker results. However, I agree that this approach is not common in other cases, there is maybe a drawback somewhere.

[simopier]

Is that the smallest you could get for this case?

[lkadz]

I tested with nl_abs_tol = 1e-8 and nl_rel_tol = 1e-6, as used in ver-1kc-1, but the testing times were excessively long, exceeding even the limit for the heavy test.

[simopier]

The legend should be put in a different spot to avoid overlap with curves (bottom right, for example).

[lkadz]

The equilibrium is much better with tighter tolerances! But I think we still have some room for improvement, especially for enclosure 2.

What was your experience playing with the tolerances and seeing their impact on the results?

Compared to the previous solving parameters, the results are noticeably better and faster. Specifically, increasing the growth factor has helped reduce the testing time, as larger time steps can be used.

In this case, I set K1=1000, but I also experimented with K1=5000. As expected, this resulted in a better chemical equilibrium in enclosure 2. However, the drawback is that the simulation takes more than 300 seconds to complete, exceeding the time limit even for the heavy test.

Moreover, this is quite strange that, while both the standard and heavy tests fail here, running them locally works without issues (no mismatches in the CSV files).

[lkadz]

Hi @simopier,
I have refined the parameters by setting nl_abs_tol = 1e-8, nl_rel_tol = 1e-6, and dtmax = 1e-4, and utilized solve_type = JFNK to achieve quicker convergence. This adjustment has resulted in better chemical equilibrium in Enclosure 2. However, the trade-off is an increase in RMSPE for the sorption law as well as a rise in mass variation.