Verification of Reaction-Diffusion Model

[smpeyres]

Hey y'all, I am trying to verify my implementation of a simple reaction-diffusion model where the electrons have some constant flux at the left boundary (@lindsayad thank you for your help with this!) and then under-go first order recombination with the solvent. I have an analytical solution for this case, and I'm trying to match my solutions.

From a linear-linear plot, they're almost identical! The zapdos model underestimates the interfacial concentration by about 0.036 mM, which I would like to improve.

I also looked at it with a log-linear plot to see if I could find out why, and turns out there's this huge numerical artifact in the middle of the domain... I expected the density to level-off on the right (no flux BC), but not the camel in the middle.

Any advice on how to fix these problems? I attached my input file and visualization script.
first_order_files.zip

[cticenhour]

Have you done a spatial convergence study to see how this discrepancy changes? I also noticed that your tolerances are relatively loose....what happens if you use the MOOSE defaults?

[cticenhour]

Not sure how to do a spatial convergence study.

You need to re-run your simulation with various levels of mesh refinement, starting somewhat coarse and ending finer. Then you take the norm of the difference between the analytical solution you have and the Zapdos result that represents your error. You can find more information in a recent Zapdos paper: https://iopscience.iop.org/article/10.1088/1361-6595/acce65.

Once you plot this (element size vs. error) on a log-log scale, you ideally will see a straight line with a particular slope, dictated by your variable finite element basis functions, showing that with increased refinement you achieve better convergence toward the exact solution. If it deviates, then it would suggest some sort of implementation error.

I used the same tolerances as the "water_only.i" test. Using defaults fails to converge, goes to dtmin.

Roughly, which combinations of values have you tried? I would first try leaving the absolute tolerance as the default and then adjusting the relative to see what happens. Tightening one or both of these tolerances might improve your result.

[smpeyres]

1. Running with 2000 elements instead of 1000 is worse -> huge dip in em value near left boundary. 500 elements about the same. 100 elements also features huge dip near left boundary. Might be implementation?

2. I've tried a couple pairs between between default and whats currently set. Using default absolute tolerance and current relative tolerance doesn't converge. Converges with 1e-7 for both at 1000 elements, but huge dip in em near left boundary (-1100 from -20).

[cticenhour]

Puzzling behavior, and it looks to me like your solution isn't "properly" converged. You're working with relatively small linear concentrations in a very small domain, but effectively not using the manual Zapdos scaling (set to 1.0), so perhaps setting automatic_scaling=true in the Executioner block in your input will help improve convergence with tightened criteria and get more consistent results.

I will not have time today to actually run this (I've only been looking at the text files thus far), but I should have time over the weekend and can give a more informed perspective.

[smpeyres]

Automatic scaling didn't converge with the loose tolerances I had originally nor with the tolerances set to 1e-7. Technically converges if I switch pjfnk to newton, but still huge dip near left boundary, None of the available integration schemes helped for newton and pjfnk.

[smpeyres]

I tried the manual Zapdos scaling of 1e-6 and changed my xmax to 1, and the solution looks much better! The flux isn't right.

Not sure how to make the units for my neumann work now (half the reason I didn't use the manual scaling at first...). The 0.0104 value was supposed to be mol/m2*s. I tried dividing the value by 1e-6 to hopefully match the scaling, and it matches the analytical result pretty well except for some bumps near the right boundary and the interfacial concentration being under-estimated by 0.0171 mM, or 1.91%.

Here is the executioner I used:

[Executioner]
type = Transient
end_time = 10
steady_state_detection = true
automatic_scaling = true
petsc_options = '-snes_converged_reason -snes_linesearch_monitor'
solve_type = newton
petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_shift_amount'
petsc_options_value = 'lu NONZERO 1.e-10'
nl_rel_tol = 1e-7
nl_abs_tol = 1e-7
dtmin = 1e-15
l_max_its = 20
scheme = crank-nicolson
[TimeStepper]
type = IterationAdaptiveDT
cutback_factor = 0.4
dt = 1e-11
growth_factor = 1.2
optimal_iterations = 30
[]
[]

Bumpiness gets much worse if my tolerances are 1e-8.

[csdechant]

Good morning. Based on your input file, you are trying to verify a Helmholtz equation (a diffusion problem with a first-order reaction source term). It looks like you have a time derivative for the electron, but your analytical solution doesn't have a time dependence?

If I remember correctly, the Helmholtz equation is ill-posed without pinning the equation. If it is time-dependent, then you need the correct initial condition at t=0, and if is time-independent, then you need to pin the solution at the walls with non-trivial derivatives and/or constant values.

What is the equation you are using for verification? Also, do you still have these convergence issues when you substitute the Zero Neumann BC with a Dirichlet BC?

[smpeyres]

I am doing an time-dependent problem to take advantage of the CRANE description for the reactions; I want to eventually expand the network and CRANE is easier for this. The analytical solution is a solution for the steady-state profile (an exponential) decay.

Pinning the equation will be hard without over-prescribing I feel like. The only idea I have to is to pin the right boundary to some very low value. Don't want to pin anything else. I also don't have a useful initial condition to help.

Let me try a very low value at the right edge. I tried -25 and -30 as a DircihletBC on the right but neither look great. Do I need another BC to say zero flux as well? Adding an extra NeumannBC of zero at the right doesn't help.

[csdechant]

For the analytical solution, I was more asking how was it derived. Based on the input file and python script, I see the following:

For the equation $\frac{\partial n_{e}}{\partial t} - \nabla\cdot (D \nabla n_e) + Kn_e =0$, where $D = 3e-10 \frac{m^{2}}{s}$ and $K = 4.5e5 \frac{1}{s}$, there is the unique solution of $n_e = \frac{G}{\sqrt{K*D}}\text{EXP}(-\sqrt{K/D} * x)$ for the conditions of $\nabla n_e \cdot \textbf{n}=G$ at $x = 0$ and $\nabla n_e \cdot \textbf{n}=0$ at $x = 1e\text{-}6m$.

My question is is that statement true (in particular, is that the correct equation and conditions used to derive the analytic solution and is there a unique solution if only Neumann BCs are used)?

For the value for the DircihletBC, it should be based on the analytic solution (what is the analytical solution at $x = 1e\text{-}6m$?). Just keep in-mind to convert the analytical solution to logarithmic form when inputting it as a DircihletBC.

[smpeyres]

I see what you mean... No, the equation was derived assuming the flux was zero at infinity.

I will try a variation where it is a steady state solve using the reaction kernels rather than CRANE. Speaking of which, I should verify if CRANE is strictly a transient code... if so, it probably could be made more flexible.

I am a little nervous about "over-prescribing" values, especially when I build complexity.

[csdechant]

Ok, does the Zapdos results start to better match the analytical solution for zero flux condition as you increase the domain (what is the $L_{2}$ error for $x$(end) = $1e\text{-}6m$, $1e\text{-}3m$, $1m$, for example). Just keep in mind that since Zapdos solves the densities in logarithmic from, so the solve can be unstable as the density approaches zero.

I will try a variation where it is a steady state solve using the reaction kernels rather than CRANE. Speaking of which, I should verify if CRANE is strictly a transient code... if so, it probably could be made more flexible.

Sorry, this is for my own information. I am a little confused by this statement. CRANE's Reactions Action Block supplies the reaction kernels (constant or energy dependent first-, second-, and third-body reactions) that are in CRANE in a user friendly way. Did CRANE's Reactions Action recently add some type of time dependent kernels? Can users no longer have steady-state reaction networks with CRANE?

[smpeyres]

I was under the impression (not sure by who) that CRANE could only do transient problems. I could be very wrong with this assumption; I've never seen a Steady problem use the CRANE Reactions block.