Calling discharge() on Well results in a ZeroDivisionError

[Huite]

import ttim

# %%

ttim_model = ttim.ModelMaq(
    kaq=[5.0],
    z=[10.0, -10.0],
    c=[],
    topboundary='conf',
    Saq=[0.001],
    Sll=[],
    phreatictop=True,
    tmin=0.01,
    tstart=0.0,
    M=10,
)
ttim_well_0 = ttim.Well(
    model=ttim_model,
    xw=0.0,
    yw=0.0,
    tsandQ=[(1.0, 500.0), (5.0, 0.0)],
    rw=0.1,
    res=0.0,
    layers=0,
    label=None,
    rc=0.0,
    wbstype='slug',
)
ttim_model.solve()

observation_peilbuizen_0=ttim_model.head(
    x=0.5,
    y=0.5,
    t=[1.5, 3.5, 5.5],
)

This runs fine, but running ttim_well_0.discharge(t=1.01) results in:

---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
[c:\tmp\timtests\ttim_discharge.py](file:///C:/tmp/timtests/ttim_discharge.py) in line 35
     [30](file:///c%3A/tmp/timtests/ttim_discharge.py?line=29) ttim_model.solve()
     [32](file:///c%3A/tmp/timtests/ttim_discharge.py?line=31) observation_peilbuizen_0=ttim_model.head(
     [33](file:///c%3A/tmp/timtests/ttim_discharge.py?line=32)     x=0.5,
     [34](file:///c%3A/tmp/timtests/ttim_discharge.py?line=33)     y=0.5,
     [35](file:///c%3A/tmp/timtests/ttim_discharge.py?line=34)     t=[1.5, 3.5, 5.5],
     [36](file:///c%3A/tmp/timtests/ttim_discharge.py?line=35) )
---> [38](file:///c%3A/tmp/timtests/ttim_discharge.py?line=37) ttim_put_0.discharge(t=1.01)

File [c:\src\ttim\ttim\element.py:203](file:///C:/src/ttim/ttim/element.py:203), in Element.discharge(self, t, derivative)
    [201](file:///c%3A/src/ttim/ttim/element.py?line=200)     s = np.sum(s[:, np.newaxis, :, :] * self.aq.eigvec, 2)
    [202](file:///c%3A/src/ttim/ttim/element.py?line=201)     s = s[:, self.layers, :] * self.model.p ** derivative
--> [203](file:///c%3A/src/ttim/ttim/element.py?line=202)     rv = invlapcomp(time, s, self.model.npint, self.model.M, 
    [204](file:///c%3A/src/ttim/ttim/element.py?line=203)                     self.model.tintervals, self.model.enumber, 
    [205](file:///c%3A/src/ttim/ttim/element.py?line=204)                     self.model.etstart, self.model.ebc, self.nlayers)
    [206](file:///c%3A/src/ttim/ttim/element.py?line=205) return rv

File [~\.conda\envs\main\lib\site-packages\numba\cpython\numbers.py:1068](https://file+.vscode-resource.vscode-cdn.net/c%3A/tmp/timtests/~/.conda/envs/main/lib/site-packages/numba/cpython/numbers.py:1068), in complex_div()
   [1066](file:///c%3A/Users/bootsma/.conda/envs/main/lib/site-packages/numba/cpython/numbers.py?line=1065) bimag = b.imag
   [1067](file:///c%3A/Users/bootsma/.conda/envs/main/lib/site-packages/numba/cpython/numbers.py?line=1066) if not breal and not bimag:
-> [1068](file:///c%3A/Users/bootsma/.conda/envs/main/lib/site-packages/numba/cpython/numbers.py?line=1067)     raise ZeroDivisionError("complex division by zero")
   [1069](file:///c%3A/Users/bootsma/.conda/envs/main/lib/site-packages/numba/cpython/numbers.py?line=1068) if abs(breal) >= abs(bimag):
   [1070](file:///c%3A/Users/bootsma/.conda/envs/main/lib/site-packages/numba/cpython/numbers.py?line=1069)     # Divide tops and bottom by b.real
   [1071](file:///c%3A/Users/bootsma/.conda/envs/main/lib/site-packages/numba/cpython/numbers.py?line=1070)     if not breal:

ZeroDivisionError: complex division by zero

[mbakker7]

I have two questions:

1. What did you have in mind with a slug test where a slug of 500 is extracted at t=1 and then another slug of 0 is extracted at t=5. What do you mean with a slug of 0?
2. A slug is extracted but the caisson has radius 0 (rc=0). That is not possible I think. It could very well be that that is where the divide by zero comes from. Please note that when rc>0 there is no error (although point 1 remains and I am not quite sure what you are trying to do).

[mbakker7]

Two suggested changes:

1. Clarify the docs that Q is the volume taken instantaneously out of the well when wbstype='slug'
2. Make sure rc > 0 when wbstype='slug'

[Huite]

I must admit I hadn't looked very carefully at the distinction between "slug" and "pumping" -- I now also realize I had set the wrong default value in the QGIS plugin.

What I wanted to compute here was pumping, from t=1.0 to t=5.0; a slug of 0 is indeed meaningless. I could've figured this out from reading the existing docs better. Having said that, being forced to think about it for a moment, it makes sense that the caisson radius should be a non-zero number. The Well docstring currently does mention everything, but it's extremely terse.

What would likely help a reader is to add some prose in the general description about the mechanism, not just in the argument list. E.g. what is potentially somewhat confusing is that the discharge extracted may be "in conflict" with the well radius, since the discharge could be far larger than the wellbore volume. Clarifying how both the radius and the caisson radius are used in the actual computation could clear that up.