Add reservoir state variable

[AnneMortensen]

Okay, so I have started on this now, and have a couple of questions.

1. While experimenting, I have gone back to my old "Annes Playground" (Offline copy, I think) where I am not using our shared "mod_funcs". Is there a better way to do this, or should I change directly in mod_funcs? (Im a little afraid to do so, I must admit)

2. Is there a good way to get the time difference (tt/nt) without having to have them as arguments for "rates"? Or is "rates" already being differentiated?

My thoughts are a mass balance like this: (the v with a line above it is the volumetric flow rate)

Where deta t is the time difference between two subsequent concentrations, and each concentration (cin=clin in model) is used for the next calculation. I have tried this, but it does not work properly (larger reservoir gives larger liquid phase concentration):

v_res>0 makes sure I dont divide by 0, c_res = clin is meant to be the large arrow shown in the first picture, and then the calculation of the new clin at the last line.

I have uploaded my files it under "Applications/modelling reservoir TEST", because this became far too complicated to understand without the actual file (still might be hard to understand what I mean even with the file). I know that multiple versions of the model is undesired, but I also did not want to push changes that do not work yet.

So, am I as far off as I feel that I am?

[sashahafner]

1. While experimenting, I have gone back to my old "Annes Playground" (Offline copy, I think) where I am not using our shared "mod_funcs". Is there a better way to do this, or should I change directly in mod_funcs? (Im a little afraid to do so, I must admit)

A new branch would be a good way to work on this. Then once it is sorted out you can merge it back with main.

Or, just work right in main. Don't worry about trying things that do not work because we can always under commits.

[sashahafner]

I have uploaded my files it under "Applications/modelling reservoir TEST", because this became far too complicated to understand without the actual file (still might be hard to understand what I mean even with the file). I know that multiple versions of the model is undesired, but I also did not want to push changes that do not work yet.

Go ahead and make the changes in the tric-fil-mod/mod_funcs.py and push them. Otherwise it is too hard to see what has changed.

I will add some suggestions for how to proceed next.

[sashahafner]

I think in general you need changes in three places:

1. rates function needs to return derivative for new state variable
2. initial state variable array needs to include one more value for reservoir mass variable
3. extract reservoir compound mass from out

These are around here:

1. https://github.com/AU-BCE-EE/tric-fil-mod/blob/6734bf4fc2e1c082eac4272046f43919c20867a2/mod_funcs.py#L129
2. https://github.com/AU-BCE-EE/tric-fil-mod/blob/6734bf4fc2e1c082eac4272046f43919c20867a2/mod_funcs.py#L217
3. https://github.com/AU-BCE-EE/tric-fil-mod/blob/6734bf4fc2e1c082eac4272046f43919c20867a2/mod_funcs.py#L231

There will be some more important details, like sorting out the exact equation for the derivative, and creating the initial state variable value depending on the input arguments. I think you are on the right track with both of these but I am not completely sure.

For the derivative, a mass balance approach is correct, yes. I would base it on mass not concentration, partially because that is what is already is used in rates (or really it is probably g per m2 cross-sectional area). So something like this:

dmcr / dt =  c_l[nc - 1] * liq_flow - c_res * liq_flow = (c_l[nc - 1] - c_res) * liq_flow

I think we can assume a constant reservoir volume, and that volume will need to be passed to rates to get c_res from the state variable mcr, the mass of compound in the reservoir.

Now I'll look at your question above again.

[sashahafner]

2. Is there a good way to get the time difference (tt/nt) without having to have them as arguments for "rates"? Or is "rates" already being differentiated?

Hmm, I don't see we would need this in rates. The ODE solver function solve_ivp does all the integration. Our rates function just has to return derivatives.

[AnneMortensen]

Add if statement to make v_res = 0 possible https://github.com/AU-BCE-EE/tric-fil-mod/blob/dbfdac31c15273c994b43dd288edf4022e8670ae/mod_funcs.py#L95C7-L95C25

[AnneMortensen]

v_res needs to be m3 pr m2 cross sectional area of the filter.

[sashahafner]

@sashahafner check L137, can't we put 3 derivatives together with concatenate?

[sashahafner]

Also 229

[sashahafner]

Do we need input reservoir initial concentration?

[sashahafner]

@sashahafner check L137, can't we put 3 derivatives together with concatenate?

Done now. Still a bit clunky, need to convert to array before concatenating.

[AnneMortensen]

hm, not sure it is working though. I am trying to run a script that worked before. So I will refer to an old version of mod_funcs until further notice.

[sashahafner]

Will take another look. I thought it worked in a test.

[sashahafner]

Still seems OK to me. Can you run demo 16? Runs for me with no problems. And with breakpoint() at L 236 before the concatenate call to check, I see

(Pdb) mcg.ndim
1
(Pdb) mcl.ndim
1
(Pdb) np.array([mcr]).ndim
1

All 1 dimension. Also at L 145:

(Pdb) dmg.ndim
1
(Pdb) dml.ndim
1
(Pdb) np.array([dmcr]).ndim
1

Do you know which line gives the concatenate dimension error?

[AnneMortensen]

Demo 16 only runs for me, when recirc is "False", because then it never enters the loop on line 96. If I put v_res = 0.001 and counter and recirc are true, I get a similar error to the one I got before. I think that it is line 145 that gives the error. In any case, it is the last line that prints before the error.

[AnneMortensen]

I made it so that the model ignores v_res if recirc is False, because physically it does not make sense to have a reservoir for the liquid phase if it is not recirculated. Maybe we should make it give an error rather than just ignoring v_res? https://github.com/AU-BCE-EE/tric-fil-mod/blob/3a0c2d51d932ff1e054b0f4757ee3cd5b433dead/mod_funcs.py#L96C1-L108C71