Add Multi Channel Transport Model (MCT)

[schmoelder]

This PR implements a modified 2D-GRM to model plant root systems as described here

To-do

• Define Interface (@jazib-hassan-juelich)
• Implementation
  • Take 2D-General Rate Model as template (rename and register classes)
  • Adapt sparse matrix structure
    https://github.com/modsim/CADET/blob/8fe0d0203a94242a0c51032704ab3585b14a1caf/src/libcadet/model/parts/TractorConvectionDispersionOperator.cpp#L1122
  • Add fluxes in residual
• Documentation (@hannahlanzrath, @lieres)
  • Model
    • Text + Equationss
    • Update figure
  • Interface
• Add unittest
  • add (dense) AD to enable AD Jacobian tests
• Check if other changes in 2D-GRM have been made since creation of this branch (e.g. section dependent velocities)
• Update exchange areas (for non-concentric channels)
• Fix build for Ubuntu and MacOS

[sleweke-bayer]

Where does the name Tractor come from? Maybe a more neutral name that can be understood by non-insiders would be fitting?
Suggestion: CompartmentTransportModel

[Immudzen]

I agree with giving it a descriptive name. That also makes it easier to use the name or if it appears in a presentation, paper etc.

[schmoelder]

Oh yes, it's just a working title because @hannahlanzrath is still checking with her supervisor what to call it.

[schmoelder]

Do we need separate function for it or can we just read the parameter like any other model parameter?

[jazib-hassan-juelich]

I don't think we will have the decay terms here. It's easier to use the existing reactions interface.

agreed, but I just added the general framework and as soon as I have information from Hannah about the actual interface then I will made modifications accordingly.

[jazib-hassan-juelich]

Do we need separate function for it or can we just read the parameter like any other model parameter?

In my view, adding a specific configuration function to read the parameters of this particular model makes the tracibility of the current implementation a bit straight forward. Just my opinion. Else we can read them like the other parameters.

[sleweke-bayer]

That will work as a first shot but is going to be inefficient.
I'd prefer setting the sparsity pattern based on the actual exchange matrix structure (i.e., only add non-zeros in the pattern for compartments that actually communicate).
To this end, setSparsityPattern() can be called in configure() instead of configureModelDiscretization().

[schmoelder]

I'd prefer setting the sparsity pattern based on the actual exchange matrix structure (i.e., only add non-zeros in the pattern for compartments that actually communicate).

We fixed that by only adding the pattern only if the entry in the matrix is non-zero.

To this end, setSparsityPattern() can be called in configure() instead of configureModelDiscretization().

@jazib-hassan-juelich @jayghoshter yesterday, you had a look at the order of calling the various configure functions. Any comments on that?

Other than that, we're still awaiting tests from @hannahlanzrath but I hope it works now.

We also implemented the exchange matrix s.t. it also would work for multiple components.

assuming e_{i,j,k} describes the transport of component k from compartment i to j, the matrix has the following form (assuming a 2 component system):
[
0, 0, e_{0,1,0}, e_{0,1,1}, e_{0,2,0}, e_{0,2,1}, ...
e_{1,0,0}, e_{1,0,1}, 0, 0, e_{1,2,0}, e_{1,2,1}, ...
...
]

@Immudzen also suggested we could take a similar approach to how we setup the connections, meaning that we only specify [compartment_from, compartment_to, component_index, exchange_rate].
This would simplify the user interface but I suggest we leave that for the python side to figure out.

[schmoelder]

To this end, setSparsityPattern() can be called in configure() instead of configureModelDiscretization().

@jazib-hassan-juelich @jayghoshter yesterday, you had a look at the order of calling the various configure functions. Any comments on that?

ping @jayghoshter @jazib-hassan-juelich

[jazib-hassan-juelich]

To this end, setSparsityPattern() can be called in configure() instead of configureModelDiscretization().

@jazib-hassan-juelich @jayghoshter yesterday, you had a look at the order of calling the various configure functions. Any comments on that?

If I remember correctly, we first tried to implement the exchange matrix in the configure function, but then we faced the problem in the residualImpl Function in Line that the exchange matrix had no entries which caused the simulation to fail. Due to this issue we defined the exchange matrix configureModelDiscretization() so that when this function is called exhange matrix is already defined.
Alongwith @Immudzen and @jayghoshter I looked deep into this issue and it relates to the order of calling these functions and when we try to change the order then broke the code. @jayghoshter can add more to this.

[jayghoshter]

Not much more to comment. As @jazib-hassan-juelich wrote, configureModelDiscretization() is called before configure() in ModelBuilderImpl.cpp:183.

if (!model->configureModelDiscretization(paramProvider, *this) || !model->configure(paramProvider))

I think the issue was that setSparsityPattern() is called in configureModelDiscretization and needs _exchangeMatrix to be defined.

[sleweke]

This will always be active. In fact, it doesn't really matter. So no need for this function to be a template. Just choose between double and active.

[sleweke]

I think the issue was that setSparsityPattern() is called in configureModelDiscretization and needs _exchangeMatrix to be defined.

Well, then don't call it at this point 😛 It makes more sense to read the exchange matrix in configure(). Instead of calling setSparsityPattern() in configureModelDiscretization(), just call it in configure() when everything is known.

[sleweke]

We have changes in two different residuals (is this correct? Could already be accounted for in _exchangeMatrix). In each change, only one state vector item is involved. Hence, we should see two changes in the Jacobian:

• d resCur_orig / d yCur_orig = exchange_orig_dest_comp
• d resCur_dest / d yCur_orig = -exchange_orig_dest_comp

We also need to take into account that the second argument of centered() refers to the offset on the main diagonal of the matrix. That is, centered(x,y) refers to matrix item mat_{x, x+y}.

This leads to the code

_jacC.centered(offsetCur_orig, 0) += static_cast<double>(exchange_orig_dest_comp);
_jacC.centered(offsetCur_dest, static_cast<int>(offsetCur_orig) - static_cast<int>(offsetCur_dest)) -= static_cast<double>(exchange_orig_dest_comp);

Since resCur_orig corresponds to yCur_orig, we get diagonal 0. For the second, we need the offset from resCur_dest to yCur_orig. From our current position (diagonal 0, equivalent to resCur_dest), we go back to the origin (-offsetCur_dest) and from there to yCur_orig (+offsetCur_orig). Note that we need to convert to int because unsigned int does not allow for negative offsets.

[sleweke]

General remarks:

• Please use consistent indentation (tabs?).
• What about some unit tests (e.g., analytic Jacobain vs. finite differences)?

[schmoelder]

sorry, I made a mess while trying to rebase. I will clean up and force push.

[sleweke-bayer]

@hannahlanzrath has put some test cases on sciebo(?). If those work correctly, we can rebase and merge.
However, adding some unit tests before merging would be great.

[schmoelder]

We need to make sure to also implement the fix for dynamic flow rates provided in #101.

[schmoelder]

Hi Hannah,

sorry for the long silence.

If we look at the to-dos, we're only missing a unit test. Maybe we can get this done next Friday during our Core-Dev Meeting. Other than that, we're good to go.

[hannahlanzrath]

MCT_demo.zip

[hannahlanzrath]

Update on the Unittest Issue

A Unittest file (test->MultiChannelTransportModel.cpp) was added, and a simple test case created, that compares the outlets of two seperate channels with the same configuration.

• From there on, more test cases can be added (I think I could do that, maybe with a little help of someone knows the language better)
• The test cases from the Jacobian must be modified so they use the MCT file (? I am not sure how that would be solved)

Also we rebased as a preparation to merge into the main branch

[hannahlanzrath]

What would be the right way to pull now to avoid merge conflicts and just get the current version?

[schmoelder]

If you don't have any local changes, try the following:

git fetch                                              # This updates the local git repo with information about the upstream branches
git reset --hard origin/tractor        # This will set the local branch to exactly the state of the remote branch.

[hannahlanzrath]

Today Jan and I worked on the Unittests:

• We added a Unittest that Compares the Output of two simulations Running Exchange from Channel 1 to Channel 2 and vice versa
• In JasonTestModels.cpp we added the case MULTI_CHANNEL_TRANSPORT with NCHANNEL and CHANNEL_CROSS_SECTION_AREAS for the Jacobian Tests and added them in test->MultiChannelTransportModel.cpp (some of those FAIL)
• We added the MCT as well for the LWE specifying NCHANNEL=1 to get the other tests working
• In test->MultiChannelTransportModel.cpp we added a Unittest for backwards flow based on the one for the other models, butit fails -> We tested forwards and backward flow and looked at the output graphically, for the backward flow it shows artifacts in the beginning even with high ncol and abs/rel tolerances for the solver. Ideas?

TODO for DevCall Nov:

• Find out why the backwards flow Test does not work
• Find out why some of the Jacobian Tests do not work (Note: For the other models there were also some tests not working)
• If these are issues with cadet core and not the tests itself: Publish MCT

[hannahlanzrath]

Jan and I found out that the Jacobian tests fail for the GRM2D as well. The fowards/backwards tests do not work for any of the models, we even tested on a branch without MCT. Thus the tests not working likely has nothing to do with the implementation of the MCT but a general problem of the cadet-core. We will check the commits if we find a state in which these tests were still working.

[hannahlanzrath]

Sam added, that for the 2DGRM Algorithmic Differentiation is not working, so

• Comment out AD tests for MCT and 2DGRM

[hannahlanzrath]

Johannes and I noticed, that the Exchange Matrix is not only nchannel x nchannel but can be defined for every component, so

• Exchange Matrix Documentation in Interface should say size nchannel x nchannel x ncomp

[jbreue16]

Ive fixed the MCT backwards flow and updated the tests:
We now have tests for simple logics: fwd vs bwd flow; fwd vs bwd exchange; two channels with same inlet yield same result.

We should add a test against a reference solution. @hannahlanzrath I think its alright to use the LRM without porosity to generate one. Maybe @sleweke can generate an analytical one using CADET-semianalytic sometime.

[ronald-jaepel]

Heads up:

Before merging, this branch would have to be rebased anyways, so we went ahead an rebased it now to get it to play nicely with the update build instructions. This will, however, break workflows for people using this branch with the old build instructions (I think). Therefore I haven't force pushed the rebased branch onto this branch but create a new one. Let's discuss in which order to merge these changes Wednesday in the office.

[schmoelder]

Congrats! 🎉

[hannahlanzrath]

thanks a lot for all your collective work on this guys, you're the coolest people 😎👌