Spurious concentrations in geochemistry

[smolins]

I am confident that PRs #290, #291, and #292 successfully address the parallel issue with transport that was affecting reactive transport simulations. (The other issue would be the overshoot issue #286 due to interpolation of saturation when subcycling that is still pending, which can be temporally addressed by not subcycling.)

Unfortunately, it does not completely resolve the problems with reactive transport simulations. Back to running the hillslope problem (Molins et al 2022 Water Resources Research) but just for a 0.1 days, capping the time step to 70 s (to avoid subcycling) and plotting every 70 seconds too, see input file attached: hillslope_calcite_crunch_sigmoid_100s_xml.txt

The issue does not seem to be related to parallel runs although parallel runs seem to provide a clue that something is wrong with the very first cell in the subdomain handled by each MPI process, as the spurious concentrations appear at regular interval across the domain. I wonder if something is done differently with the first cell (at initialization or after solving transport) that affects (not so much the total concentration tcc) but the primary ion concentration (or another variable handled via Alquimia). The issue does not appear for non-reactive tracers (1 and 2), only with reactive species.

Time step 1, single processor (issue is at top left corner)

Time step 1, parallel run (n=12) (issue at top left corner repeats itself downstream)

Time step 10, single processor

Time step 10, parallel run

[smolins]

@dasvyat: could you please have an early look at this and provide some clues as to what may be happening? Thanks.
It has me confused that it is just for the one cell at the beginning of the subdomain. The concentrations should all be the same in these examples (i.e. red).

[smolins]

data and mesh directories for this example are available in ats-demos under 13_integrated_hydro_reactive_transport

[levuvietphong]

@smolins: Are the plots you show from a parallel run? If yes, how many cores did you use?

[smolins]

Plots 1 and 3 are for single processor runs. They show the spurious concentrations on the top left corner.

Plots 2 and 4 are for parallel runs using 12 processors. They show the spurious concentrations on the top left corner and then twice more downstream. This is a zoomed-in view. If showing the entire domain, one would see this pattern 12 times in total.

[levuvietphong]

@smolins I re-run your input file, and it appears that:

• The problem may be experiencing a dry condition on the surface since surface-water_source or rainfall=0 for the test period (0.1 day)
• The colorbar range is very narrow (e.g., for Ca++, the difference between red and blue is ~1e-7), meaning the variation is small. This could be due to the dry conditions, and the tcc calculation may be sensitive to this dry condition.

I added some rainfall at the beginning, and this seems eliminate the weird pattern. Can you play with the rain and verify this?

[dasvyat]

I've re-run this case on multiple cores. There is a parallel pattern that indicates a bug versus dry cells or boundary conditions.
I'll take a look at this problem.

[dasvyat]

@smolins . Is it possible that Alquimia is doing smth strange with the first cell and some quantities are not initialized. This some additional output that I've got for this chemical engine:

cell 0
aqueous_components 6.9448e-05 0.0020636 0.000997075 0.000997075 9.97075e-11
saturation 0.279566
volume 0.5
porosity 0.25
den 1000
temp 25
mineral_volume_fraction 0.01
mineral_rate_cnst 0
mineral_specific_surface_area 1
cell 0
****Alquimia out -1.117e-05 0.00198136 0.000996266 0.000997075 9.97075e-11
mineral_volume_fraction 0.01
mineral_rate_cnst 0
mineral_specific_surface_area 1

cell 1
aqueous_components 6.9448e-05 0.0020636 0.000997075 0.000997075 9.97075e-11
saturation 0.279963
volume 0.5
porosity 0.25
den 1000
temp 25
mineral_volume_fraction 0.01
mineral_rate_cnst 0
mineral_specific_surface_area 1
cell 1
****Alquimia out 1.91651e-05 0.00201258 0.000996709 0.000997075 9.97075e-11
mineral_volume_fraction 0.01
mineral_rate_cnst 0
mineral_specific_surface_area 1

cell 2
aqueous_components 6.9448e-05 0.0020636 0.000997075 0.000997075 9.97075e-11
saturation 0.280364
volume 0.5
porosity 0.25
den 1000
temp 25
mineral_volume_fraction 0.01
mineral_rate_cnst 0
mineral_specific_surface_area 1
cell 2
****Alquimia out 3.52199e-05 0.00202895 0.000996865 0.000997075 9.97075e-11
mineral_volume_fraction 0.01
mineral_rate_cnst 0
mineral_specific_surface_area 1

[dasvyat]

In spite very similar input values provided by ATS for the first cell (cell=0) Alquimia produces negative values for the first component. @smolins does it look suspicious for you?

[smolins]

That is kind of what I suspected in the Issue description. I wonder whether the initial guesses supplied to Alquimia in the first cell are not good and the solution for whatever reason is not good (behavior then is not captured by Alquimia). Is there a way to print out other data provided into Alquimia, not just total concentrations and mineral stuff? I mean all other stuff passed back from the Alquimia initialization step must be stored and supplied back to Alquimia for time stepping, not just total concentrations.

…

On Fri, Feb 7, 2025 at 2:55 PM Daniil Svyatsky ***@***.***> wrote: In spite very similar input values provided by ATS for the first cell (cell=0) Alquimia produces negative values for the first component. @smolins <https://github.com/smolins> does it look suspicious for you? — Reply to this email directly, view it on GitHub <#293 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADXOTGRTIJBDIRUK3BOJNKD2OU2WXAVCNFSM6AAAAABWGMQKMWVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDMNBUGI3TOMBZGY> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[lipnikov]

Are you using the log formulation?

I think if you initialize cells 0 and 1 to the same values, the result will be the same. But it would be good to verify this.

If you switch to the native chemistry, you may get more useful info about beaker initialization in cell 0.

[smolins]

Crunch always uses the log formulation. For the same values going in, the same results should be obtained (except for the small difference in the saturation). But it is also possible something else happens and we're getting just garbage not an actual solution. Seeing alquimia auxiliary data (not output auxiliary) going in would help.

…

On Sun, Feb 9, 2025 at 11:36 AM Konstantin Lipnikov < ***@***.***> wrote: Are you using the log formulation? I think if you initialize cells 0 and 1 to the same values, the result will be the same. But it would be good to verify this. If you switch to the native chemistry, you may get more useful info about beaker initialization in cell 0. — Reply to this email directly, view it on GitHub <#293 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADXOTGSPIKHD742SSO2FWFL2O6U4FAVCNFSM6AAAAABWGMQKMWVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDMNBWGU2TGNBRHE> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[dasvyat]

What is convincing me that it could be Alquimia issue is that if I comment out just Alquimia update - chem_engine_->Advance - this pattern disappears. Species are processed through the whole chemistry workflow but they are not touched by chem_engine_.