Add Multi Channel Transport Model (MCT)

[schmoelder]

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

To-do

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

[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.

[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]

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]

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.

[hannahlanzrath]

After some name changes of variables to match the MCT model characteristics, the documentation has to be updated.

• [x] Remove everything related to film diffusion
• [x] Change variable names / delete unneeded variables
• [ ] Add cross-section scheme to show the differences to original model with radial zones

Update!!! RESOLVED, check next comment for current TODOS:

I am not entirely sure about how the variable names in cadet are called and could use some clarification. In src/libcadet/model/MultiChannelTransportModel.hpp lines 264-265

		virtual unsigned int numAxialCells() const CADET_NOEXCEPT { return _disc.nCol; }
		virtual unsigned int numRadialCells() const CADET_NOEXCEPT { return _disc.nRad; }

was updated to

		virtual unsigned int numPrimaryCoordinates() const CADET_NOEXCEPT { return _disc.nCol; }
		virtual unsigned int numSecondaryCoordinates() const CADET_NOEXCEPT { return _disc.nRad; }

But since in the documentation the variables are called ncol / nrad, and the return _disc.nCol;/ return _disc.nRad; is it still called with the same name from the python interface using ncol / nrad in the dictionary?

I think the names in the documentation should match the new primary / secondary coordinates nomenclature or something along these lines, just not sure if I can simply change it or if more code changes are needed.

[hannahlanzrath]

TODO:

In MultiChannelConvectionDispersionOperator.cpp:

Change Parameters Names to match Documentation

• [x] Change NRAD to NCHANNEL

(- [x] Change RADIAL_DISC_TYPE to CHANNEL_TYPE ) --> Solved by defining Cross Section of each channel separately keep options 'Equivolume' & 'User_defined'! (- [x] Remove option 'Equidistant' from CHANNEL_TYPE (former RADIAL_DISC_TYPE) and set 'Equivolume' as default )

https://github.com/modsim/CADET/blob/8ef7043dcaa6948a745b117235bd4d42228e37ae/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp#L757

Change the way the areas of the Channels is defined

CHANNEL_CROSS_SECTION_AREAS is a vector of length NCHANNEL filled with scalars that define the absolute value of the cross section area for each channel

• [x] Change RADIAL_COMPARTMENTS to CHANNEL_CROSS_SECTION_AREAS depending on the method we choose to define the cross section area of each channel:

https://github.com/modsim/CADET/blob/8ef7043dcaa6948a745b117235bd4d42228e37ae/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp#L765

Reference for COL_RADIUS:

https://github.com/modsim/CADET/blob/8ef7043dcaa6948a745b117235bd4d42228e37ae/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp#L746

• [x] Check that the Mass Balance between Channels with differing Cross Section Areas is correct, compare to 2DGRM, if needed ADD Factor that weighs in Channel Volumina/Area Differences

https://github.com/modsim/CADET/blob/8ef7043dcaa6948a745b117235bd4d42228e37ae/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp#L1051

Other

• [x] Remove COL_POROSITIES from Code, as it is not needed nor used in the MCT

https://github.com/modsim/CADET/blob/8ef7043dcaa6948a745b117235bd4d42228e37ae/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp#L747

@hannahlanzrath

• [ ] Provide another Jupyter Notebook with test cases Including mass balance test for exchange between channels with differing volumes
• [x] Update documentation, units are missing
• [ ] Documentation for external functions: "radial" coordinate represents channels and is given by the respective index (i.e. r= 0 for channel 0)

[hannahlanzrath]

Hey guys, @schmoelder @jbreue16

I have a question regarding the parameters GS_TYPE , MAX_KRYLOV, MAX_RESTARTS and SCHUR_SAFETY. Are these still relevant for the MCT? As cadet does not drop an error when I assign nonsensical values such as strings to those parameters and Eric was not sure how they would come into play in the MCT, the question arose if they can be dropped from the documentation.

Also maybe you can help me understand the parameter SENS_PARTYPE. In VELOCITY and COL_DISPERSION it was specified, that "In case of a spatially inhomogeneous setting, the SENS_PARTYPE field is used for indexing the radial cell when specifying parameter sensitivities." Is the radial cell the same as a radial zone and is it now analog for the channels, so I can just change radial cell to channel?

Thanks!

[jbreue16]

We don't need those parameters. They're required for the block decomposed linear solver which is only used for models with pores.

"Is the radial cell the same as a radial zone and is it now analog for the channels, so I can just change radial cell to channel?" - Yes. Maybe we should also rename SENS_PARTYPE accordingly to something like SENS_CHANNEL ?

[lieres]

What would SENS_CHANNEL mean? There are no particles in the MCT model.

[hannahlanzrath]

Since SENS_PARTYPE is defined in the /input/sensitivity group (https://cadet.github.io/master/interface/sensitivities.html#group-input-sensitivity) I think it should be left as it is but maybe omitted from the MCT documentation, when we are sure it is not needed. My problem is that I do not really understand what it does exactly, but I also think we might not need it in the MCT model?

[lieres]

I don't think it makes sense to leave it in the code when no particles exist in the model.

[hannahlanzrath]

I accidentally pushed all my files instead of just the documentation changes, can someone assist me in rebasing maybe? Sorry, I tried to fix it on my own but I don't want to accidentally delete something.

EDIT: Problem is solved!

[jbreue16]

SENS_PARTYPE is (re)used for parameter sensitivities of parameters that depend on the particle type or the radial position (i.e. channel for the MCT), so this should definitely stay in the code. I was just wondering if we should rename the parameter as we have a model with radial position/channels but without particles types. I am not sure what's less confusing.

Iam not an expert with git, maybe we should wait for Johannes. Or, if you want, we can have a look at this tomorrow, shouldn't be too difficult.

[hannahlanzrath]

I managed to fix my git mistake.

[schmoelder]

Hi Hannah, what's the current status on this PR? Please give us a sign when this is ready to merge. Best Johannes

[hannahlanzrath]

Hello everyone!

I created two Jupyter Notebooks that Test and Demonstrate the Features of the MCT.

https://fz-juelich.sciebo.de/s/4caow0xsf6DLv7J https://fz-juelich.sciebo.de/s/PUKeXZWw0ltOo6e

Maybe someone can download them and double check, that this is working with the current build of the MCT. If everything is fine, we can merge.

Best Hannah

[jbreue16]

Hi @hannahlanzrath ,

I built the current MCT branch and your notebooks run fine with it.

[hannahlanzrath]

Then I think we're good to go!

[hannahlanzrath]

Hello everyone, just because everyone is off and on holiday, before I am gone for a few weeks I wanted to ask what is the current status and plan with merging the MCT?

[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.