How to change time-varying boundary condition after a model is created?

[MattBrst]

I would like to change a time-varying boundary condition (in my case: tsandhstar) after the model has been created. How can this be done?

(Originally asked here, but moved to this new discussion.)

[mbakker7]

Consider the model that was suggested in Discussion mbakker7/ttim#97, which is a cross-sectional model with a river with a leaky bottom on the left and an aquifer with a semi-confining top on the right. The head in the river goes up from 0 to 1 at time t=0, as specified with tsandhstar (note that I changed the names somewhat to avoid confusion):

import numpy as np
import ttim

ml = ttim.ModelXsection(naq=1, tmin=1e-4, tmax=1e2)

z = [5, 0, -10]
kh = 10
c = 100
Saq = 1e-3
Sll = 1e-2
tsandhstar = [(0, 1)]

riv = ttim.XsectionMaq(
    ml,
    -np.inf,  # river extends to infinitiy
    0,
    z=z,
    kaq=kh,
    Saq=Saq,
    Sll=Sll,
    c=c,
    topboundary="semi",
    tsandhstar=tsandhstar,
    name="river",
)

land = ttim.XsectionMaq(
    ml,
    0,
    np.inf,  # hinterland extends to infinity
    kaq=kh,
    z=z,
    Saq=Saq,
    Sll=Sll,
    c=c,
    topboundary="semi",
    name="landpart",
)

ml.solve()

The question now is, how can the user change the head variation after the model has been created?

Every model has an attribute aq, which has a dictionary with all the inhomogeneities called inhomdict.
Here, the model is stored in the variable ml and the river was called "river" when it was created (name="river"). We want to redefine the time variation of the head in the river, which is stored in the variable tsandhstar. For example, let's try to change the river stage back to 0 after 4 days. Make sure to solve the model again to initialize the changed inhomogeneity correctly.

ml.aq.inhomdict['river'].tsandhstar = [(0, 1), (4, 0)]
ml.solve() # don't forget!

After this change, the river stage in the model varies according to the new definition. The same works on a copy of the model. Note that to make a copy of a model, you should use deepcopy as explained in Discussion mbakker7/ttim#97.

Final comment: changing the time variation of an element (like a well or line-sink) rather than an inhomogeneity cannot be done in the same way at the moment. I think this would be useful, so I opened Issue #107.

[mbakker7]

@MattBrst if you would be interested to write the above answer up as an HowTo for the documentation, including one or two graphs, that may be a nice addition.

[MattBrst]

@mbakker7: Thank you for your support so far. Is there a way to use the land shorthand in copied models instead of having to use the ml.aq.inhomdict["landpart"] expression for further actions?

I'll illustrate my question with some more code which is a combination of mbakker7/ttim#97 and mbakker7/ttim#106:

# %% 1. Setup a model and run a first calibration
# Note how easily the inhomogeneity is defined in the first calibration run.
import numpy as np
import ttim

ml = ttim.ModelXsection(naq=1, tmin=1e-4, tmax=1e2)

z = [5, 0, -10]
kh = 10
c = 100
Saq = 1e-3
Sll = 1e-2
tsandhstar = [(0, 1)]

riv = ttim.XsectionMaq(
    ml,
    -np.inf,  # river extends to infinitiy
    0,
    z=z,
    kaq=kh,
    Saq=Saq,
    Sll=Sll,
    c=c,
    topboundary="semi",
    tsandhstar=tsandhstar,
    name="river",
)

land = ttim.XsectionMaq(
    ml,
    0,
    np.inf,  # hinterland extends to infinity
    kaq=kh,
    z=z,
    Saq=Saq,
    Sll=Sll,
    c=c,
    topboundary="semi",
    name="landpart",
)

ml.solve()

cal = ttim.Calibrate(ml)

cal.series(
    "pb",
    5,
    0,
    layer=0,
    t=[0,1],
    h=[0,1],
)

cal.set_parameter(
    "kaq",
    layers=0,
    initial=10,
    pmin=1.0,
    pmax=100.0,
    inhoms=([land])
)

cal.fit_least_squares(report=True, method="trf")
print(f"RMSE: {cal.rmse():.5f}")

# %% 3. Copy the model and run a second calibration
from copy import deepcopy
ml2 = deepcopy(ml)
ml2.solve()

cal2 = ttim.Calibrate(ml2)

cal2.series(
    "pb",
    20,
    0,
    layer=0,
    t=[10, 50],
    h=[0.1, 0.2],
)

# The easy definition would be something like this, but does not work yet:
cal2.set_parameter(
    "kaq", layers=0, initial=10, pmin=1.0, pmax=100.0, inhoms=([land])
)

# This is the workaround according to mbakker7/ttim#97, but is not very elegant:
cal2.set_parameter("kaq", layers=0, initial=10, pmin=1.0, pmax=100.0, inhoms=(ml2.aq.inhomdict["landpart"]))

cal2.fit_least_squares(report=True, method="trf")

[mbakker7]

What do you not like about the ml.aq.inhomdict["landpart"] and ml2.aq.inhomdict["landpart"]?
It is consistent.
I understand that you saved ml.aq.inhomdict["landpart"] also in the variable land, but that cannot be copied. And that would also be confusing, to store an inhom in a variable and give it a name. You don't have to store it in a variable at all. Will that reduce the confusion?

You talk about a workaround, but the code works exactly as designed, so I wouldn't call it a workaround.

[MattBrst]

Thank you for your clarification concerning the use of variables like land.

Based on your reply I figured out that adding double quotes and using the name of the inhom instead of the variable solves my question for the calibration case:
cal2.set_parameter("kaq", layers=0, initial=10, pmin=1.0, pmax=100.0, inhoms=["landpart"])

[mbakker7]

Maybe the doc string wasn't clear enough on how to define the inhoms? I think inhoms can be defined either by variable where it is stored or by name (string). Maybe we should change the example notebooks to use the strings and not store the inhoms in variables to avoid confusion? (Especially since that is the preferred way to do it.)

[dbrakenhoff]

Good idea, referring to the string is the way to go :)

[mbakker7]

I'll make those changes