Coupled flow and transport PK with more complex water balance

155 views
Skip to first unread message

Sam Shaheen

unread,
May 16, 2025, 9:17:56 AMMay 16
to Amanzi-ATS Users
Hi all,
First off, thanks for all the hard work on the new transport update. I updated two of the integrated hydro reactive transport demos to run with the new coupled flow and transport PK in the latest master and things look great (it fixed the convergence/concentration issue in the hillslope demo without having to reduce the flow timestep to avoid subcycling).

I am testing it out on a more complex simulation now, where the "flow" is a water balance weak MPC comprising Richards flow, overland flow, and canopy and snow general surface balance PKs, rather than just Richards and overland flow like the demos. 

With canopy and snow PKs included, I'm getting the error message:
what(): Evaluator "canopy-throughfall_drainage_rain" @ "reactive_transport_next" cannot be created in State. Verify (1) SetEvaluator is called or (2) name exists in state->evaluators.

I am wondering if my current PK tree setup is not well-suited for the new coupled flow and transport PK/whether coupled flow and transport is compatible with canopy and snow parameters in its current iteration. I'm attaching a screenshot of my current PK tree setup.

Thanks!
Sam


PKtree.png

Phong Le

unread,
May 16, 2025, 1:24:05 PMMay 16
to Sam Shaheen, Amanzi-ATS Users

Hi Sam,

Your PK tree looks good. I'm also running reactive transport with snow and canopy PKs, and I've included the subgrid hyporheic zone PK in the tree as well (see attached). However, I haven't tested my setup on the updated master branch, so I am not sure whether the latest update might have any new issues. I'll test it soon.

Regarding the error you show, it means that the "canopy-throughfall_drainage_rain" is being evaluated at the next time step that is related to the transport PK (@ "reactive_transport_next"). I think this is problematic, as canopy throughfall has nothing to do with the transport PK. It should be evaluated at the next step of the water balance PK only.

Were you able to run this setup successfully with the previous version of ATS?

Phong

Screenshot 2025-05-16 at 1.07.18 PM.png
--
You received this message because you are subscribed to the Google Groups "Amanzi-ATS Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ats-users+...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/ats-users/0d4dc6e2-840e-4810-8b15-071d29225ce7n%40googlegroups.com.


--
P

--
You received this message because you are subscribed to the Google Groups "Amanzi-ATS Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ats-users+...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/ats-users/0d4dc6e2-840e-4810-8b15-071d29225ce7n%40googlegroups.com.


--
P

Sam Shaheen

unread,
May 16, 2025, 1:59:51 PMMay 16
to Amanzi-ATS Users
Hi Phong,
Thanks. This was directly carrying over a PK tree that I have been using successfully in previous versions of master, but updating the flow and transport from a subcycling MPC to the new coupled flow and transport.

In the demos, I got a similar error message when I forgot to update flow and transport from a subcycling MPC (but with darcy velocity instead of canopy drainage). So I get the sense this has to do with the new coupled flow and transport and which variables are passed between flow and transport. Maybe there's a better way to structure water balance in light of this?
-Sam

Phong Le

unread,
May 16, 2025, 3:06:37 PMMay 16
to Amanzi-ATS Users
Yep, the transport PK iterates through a list of variables and checks whether each one is created @reactive_transport_next before advancing the transport simulation for each timestep. The first missing variable will trigger the error and stop the simulation.
Can you try running two versions - one with subcycling and one without - to see if you get the same error? This might help us isolate the issue. It's possibly a bug, and we'll need to debug and fix it.

Phong

Sam Shaheen

unread,
May 19, 2025, 9:08:12 AMMay 19
to Amanzi-ATS Users
Hi Phong,
Yep, I'll test it out and report back with what I find.
-Sam

Coon, Ethan

unread,
May 19, 2025, 12:18:29 PMMay 19
to Sam Shaheen, Amanzi-ATS Users

If the problem is what I suspect, it’s not exactly a bug in the code, but may be something we can deal with in a code change.  Long email incoming…

 

In your problem, I suspect you’re defining a source for transport that is dependent upon a source for flow.  Is the transport source given by “canopy-drainage” directly?  Probably yes – you’re probably doing something like “surface-water_source” = “canopy-drainage” – “surface-evaporation”, but want to prescribe a concentration of the incoming water only?  I’m going to assume this is the physics you’re solving, but correct me if I’m wrong and we can iterate.

 

The new “time tags” used in subcycling make things a bit more formal in terms of how you define temporal interpolation within a subcycled timestep.  If flow advances without subcycling (from CURRENT NEXT) and transport is subcycling, then transport needs to know all values at some intermediate time “reactive_transport_next” which is partway between CURRENT and NEXT as it subcycles between the two “outer” tags.

 

When transport is subcycled relative to flow with a source that depends upon the flow source (e.g. Q_C = C * Q_water), this means we need to define an interpolant for the flow source.  Since we typically treat the flow source as a “backward Euler/implicit” term in flows time integration, this implies that the source is “piecewise constant” using the value from the “NEXT” time.  So it’s as simple as telling the code that “canopy-drainage@reactive_transport_next” = “canopy-drainage@next”.

 

This is done through an “alias” evaluator, which effectively just says that the value at the time “reactive_transport_next” will always be equal to the value at “next”.  You can do this manually by adding an evaluator to your input file,

 

State->evaluators->”canopy-drainage@reactive_transport_next”

  evaluator type: alias

  target: canopy-drainage@next

 

Note that there are different possibilities here.  You could use a “Crank-Nicholson” time discretization for your source term in the flow equation, which implies a piecewise linear source.  In that case, you’d have to use a “temporal interpolation” evaluator to linearly interpolate your subcycled time between CURRENT and NEXT times in order to keep transport and flow consistent.

 

So the short-term fix is probably just do the above manually in your input file.

 

In the case of “water_flux” or Darcy flux (which like you discovered, needs the same thing) we now use the MPC that couples flow and transport to actually write the input spec in the code for the user – in this case, the MPC knows that the flow PK always uses backward Euler flow time integration, and so the flux is always going be backward Euler.  See

 

$ATS_SRC_DIR/src/pks/mpc/mpc_flow_transport.cc:69

 

To see where the MPC writes this into the state->evaluators list.  (Note, this is different from how it used to be, which hard-coded this as NEXT.)

 

So the question becomes, can we figure out how to be smart enough to write this evaluator into the input spec for the user?

 

To do this, the MPC coupling flow and transport would have to:

 

  1. Parse the transport input spec to figure out what flow sources are used (e.g. figure out that the source in question is “canopy-drainage”)
  2. Check the flow input spec to figure out that, yes, this in fact a part of flow’s sources
  3. Check whether flow’s source will be integrated using forward Euler (piecewise constant with value = CURRENT), backward Euler (piecewise constant with value = NEXT), or CN (piecewise linear).
  4. Write the resulting evaluator into the list.

 

2 in particular is hard, because “canopy-drainage” isn’t the name of flow’s source (that is surface-water_source).  But it is a dependency of the source.  So it might be doable.  We can add a ticket for this and Phong and I can think about it and see if we can get something to do this automatically.

 

Ethan

 

 

 

 

 

 

 

From: ats-...@googlegroups.com <ats-...@googlegroups.com> on behalf of Sam Shaheen <samuelw...@gmail.com>
Date: Monday, May 19, 2025 at 7:08 AM
To: Amanzi-ATS Users <ats-...@googlegroups.com>
Subject: [EXTERNAL] Re: Coupled flow and transport PK with more complex water balance

Hi Phong, Yep, I'll test it out and report back with what I find. -Sam On Friday, May 16, 2025 at 3: 06: 37 PM UTC-4 levuvi. . . @ gmail. com wrote: Yep, the transport PK iterates through a list of variables and checks whether each one is created @reactive_transport_next

Hi Phong,

Yep, I'll test it out and report back with what I find.

-Sam

 

On Friday, May 16, 2025 at 3:06:37 PM UTC-4 levuvi...@gmail.com wrote:

Yep, the transport PK iterates through a list of variables and checks whether each one is created @reactive_transport_next before advancing the transport simulation for each timestep. The first missing variable will trigger the error and stop the simulation.

Can you try running two versions - one with subcycling and one without - to see if you get the same error? This might help us isolate the issue. It's possibly a bug, and we'll need to debug and fix it.

 

Phong

On Friday, May 16, 2025 at 1:59:51 PM UTC-4 samuelw...@gmail.com wrote:

Hi Phong,

Thanks. This was directly carrying over a PK tree that I have been using successfully in previous versions of master, but updating the flow and transport from a subcycling MPC to the new coupled flow and transport.

 

In the demos, I got a similar error message when I forgot to update flow and transport from a subcycling MPC (but with darcy velocity instead of canopy drainage). So I get the sense this has to do with the new coupled flow and transport and which variables are passed between flow and transport. Maybe there's a better way to structure water balance in light of this?

-Sam

 

On Friday, May 16, 2025 at 1:24:05 PM UTC-4 levuvi...@gmail.com wrote:

 

Hi Sam,

 

Your PK tree looks good. I'm also running reactive transport with snow and canopy PKs, and I've included the subgrid hyporheic zone PK in the tree as well (see attached). However, I haven't tested my setup on the updated master branch, so I am not sure whether the latest update might have any new issues. I'll test it soon.

 

Regarding the error you show, it means that the "canopy-throughfall_drainage_rain" is being evaluated at the next time step that is related to the transport PK (@ "reactive_transport_next"). I think this is problematic, as canopy throughfall has nothing to do with the transport PK. It should be evaluated at the next step of the water balance PK only.

 

Were you able to run this setup successfully with the previous version of ATS?

 

Phong

 

Sam Shaheen

unread,
May 19, 2025, 8:48:55 PMMay 19
to Amanzi-ATS Users
Thanks Ethan!
That makes sense and sounds manageable. That said, I found inserting:
      <ParameterList name="canopy-drainage@reactive_transport_next">
        <Parameter name="evaluator type" type="string" value="alias"/>
        <Parameter name="target" type="string" value="canopy-drainage@next"/>
      </ParameterList>

gives me a similar error message that the canopy-drainage@next evaluator cannot be created in state. If I alias that to just canopy-drainage it will run (after doing the same for snow-depth, surface-radiation_balance, and saturation-gas), but it looks like I broke transport doing that as concentrations are 0 for all components across the domain throughout the run (I wasn't sure removing the @next was valid, just curious if it would eventually run).

Thanks,
-Sam

Coon, Ethan

unread,
May 19, 2025, 11:29:56 PMMay 19
to Sam Shaheen, Amanzi-ATS Users

Right, sorry, I gave you the wrong name.  The Tags::NEXT tag is actually the default tag, and so is the empty string.

 

So it should probably be:

 

        <Parameter name="target" type="string" value="canopy-drainage@"/>

Just to make it clearer that this is not a variable name without a tag, but a variable name at the default tag.

 

I suspect what you did,

 

        <Parameter name="target" type="string" value="canopy-drainage"/>

also does the right thing.

 

So the next question is why are you getting 0 concentrations.  I’m guessing you either 1, need to run this through the input converter (tools/input_converters/xml-1.5-master.py), and/or maybe 2, we don’t have the input converter quite right yet.  Likely the problem is your source term should provide “mole ratio” instead of “total component concentration”.

 

Reply with the xml for your (transport) source term, and I can compare to what I think should be there for master.

 

Thanks for helping us debug this – it is very new code.

 

Ethan

 

 

 

 

 

From: ats-...@googlegroups.com <ats-...@googlegroups.com> on behalf of Sam Shaheen <samuelw...@gmail.com>
Date: Monday, May 19, 2025 at 6:49 

PM
To: Amanzi-ATS Users <ats-...@googlegroups.com>
Subject: Re: [EXTERNAL] Re: Coupled flow and transport PK with more complex water balance

Thanks Ethan! That makes sense and sounds manageable. That said, I found inserting: <ParameterList name="canopy-drainage@ reactive_transport_next"> <Parameter name="evaluator type" type="string" value="alias"/> <Parameter name="target"

Thanks Ethan!

That makes sense and sounds manageable. That said, I found inserting:

      <ParameterList name="canopy-drainage@reactive_transport_next">
        <Parameter name="evaluator type" type="string" value="alias"/>
        <Parameter name="target" type="string" value="canopy-drainage@next"/>
      </ParameterList>

 

gives me a similar error message that the canopy-drainage@next evaluator cannot be created in state. If I alias that to just canopy-drainage it will run (after doing the same for snow-depth, surface-radiation_balance, and saturation-gas), but it looks like I broke transport doing that as concentrations are 0 for all components across the domain throughout the run (I wasn't sure removing the @next was valid, just curious if it would eventually run).

 

Thanks,

-Sam

On Monday, May 19, 2025 at 12:18:29 PM UTC-4 Coon, Ethan wrote:

If the problem is what I suspect, it’s not exactly a bug in the code, but may be something we can deal with in a code change.  Long email incoming…

 

In your problem, I suspect you’re defining a source for transport that is dependent upon a source for flow.  Is the transport source given by “canopy-drainage” directly?  Probably yes – you’re probably doing something like “surface-water_source” = “canopy-drainage” – “surface-evaporation”, but want to prescribe a concentration of the incoming water only?  I’m going to assume this is the physics you’re solving, but correct me if I’m wrong and we can iterate.

 

The new “time tags” used in subcycling make things a bit more formal in terms of how you define temporal interpolation within a subcycled timestep.  If flow advances without subcycling (from CURRENT à NEXT) and transport is subcycling, then transport needs to know all values at some intermediate time “reactive_transport_next” which is partway between CURRENT and NEXT as it subcycles between the two “outer” tags.

Error! Filename not specified.

Sam Shaheen

unread,
May 27, 2025, 2:44:39 PMMay 27
to Amanzi-ATS Users
Hi Ethan,
Re: 0 concentration, I ran my scripts through the input converter before starting out. That did still need the primary variable key suffix in the transport PKs to be changed from total component concentration to molar ratio in order to run successfully, so it's possible that there's other key changes that I (and the script) overlooked. That said, the two reactive transport demos I tested (hillslope_calcite_crunch_sigmoid and integrated_hydro_reactive_transport-superslab) ran quite well after I ran xml-1.5-master.py and made this change so I'm guessing the root cause of the issue is something related to a more complex water balance.

I'm attaching the transport portion of my input script-I played around a bit with molar_ratio vs. total_component_concentration in the source terms but it didn't fix the 0 concentration issue. If I have time this week, I'll play around with adding some complexity on the snow/canopy side to the demos to see if they still run well.

Thanks,
Sam
transport.xml

Ethan Coon

unread,
May 27, 2025, 3:15:07 PMMay 27
to Sam Shaheen, Amanzi-ATS Users
Yes, you're seeing the same problem Phong had, and I need to work on the input converter.  The problem is here:

      <ParameterList name="boundary conditions">
        <ParameterList name="concentration">
          <ParameterList name="BC coupling">
            <Parameter name="spatial distribution method" type="string" value="domain coupling" />
            <Parameter name="submodel" type="string" value="conserved quantity" />
            <Parameter name="regions" type="Array(string)" value="{surface}" />
            <ParameterList name="fields">
              <Parameter name="conserved_quantity_key" type="string" value="surface-total_component_quantity" />
              <Parameter name="field_out_key" type="string" value="surface-molar_ratio" />
            </ParameterList>
          </ParameterList>
        </ParameterList>
      </ParameterList>

You shouldn't write this list -- ATS's MPC for coupled transport will write it for you (if it doesn't exist).  Remove this completely, leaving an empty "boundary conditions" list.

Also, this list has changed, and is also written by the MPC, so should be removed:

      <ParameterList name="source terms">
        <ParameterList name="concentration">
          <ParameterList name="coupling">
            <Parameter name="regions" type="Array(string)" value="{surface domain}" />
            <Parameter name="spatial distribution method" type="string" value="domain coupling" />
            <Parameter name="submodel" type="string" value="rate" />
            <ParameterList name="fields">
              <Parameter name="flux_key" type="string" value="mass_flux" />
              <Parameter name="copy_flux_key" type="string" value="next_timestep" />
              <Parameter name="field_out_key" type="string" value="molar_ratio" />
              <Parameter name="field_in_key" type="string" value="surface-molar_ratio" />
            </ParameterList>
          </ParameterList>
        </ParameterList>



Try comparing your xml to, for instance, ats-regression-tests/06_transport/column_evaporation.xml which also solves integrated flow coupled to integrated transport.

Ethan







To view this discussion visit https://groups.google.com/d/msgid/ats-users/4d7ca69c-d726-4129-8e7b-be2c5965ce8an%40googlegroups.com.


--

-------------------------------------------------------------------
Ethan Coon
917-969-6831
https://www.ornl.gov/staff-profile/ethan-t-coon
-------------------------------------------------------------------

Sam Shaheen

unread,
Jun 4, 2025, 9:15:33 PMJun 4
to Amanzi-ATS Users
Hi Phong and Ethan,
After playing around with a lot of different PK and evaluator set-ups in some updated reactive transport demos, I'm confident that the issue I'm having isn't related to a more complex water balance PK or issues in the input script with the new master. Rather, it seems to be the parallelization issue that Yi reported in a separate thread. Running the simulation on 1 processor fixes the 0 concentration issue, and things generally look great beyond that once the handful of @reactive_transport_next evaluators are inserted.
Best,
Sam

Coon, Ethan

unread,
Jun 5, 2025, 9:43:04 AMJun 5
to Sam Shaheen, Amanzi-ATS Users

Ok, seems like we need to do some more parallel testing.  Thanks for checking,


Ethan

 

Reply all
Reply to author
Forward
0 new messages