diff --git a/docs/steady/02examples/circular_area_sink.ipynb b/docs/steady/02examples/circular_area_sink.ipynb index 6d99ecc..f1be9ea 100644 --- a/docs/steady/02examples/circular_area_sink.ipynb +++ b/docs/steady/02examples/circular_area_sink.ipynb @@ -52,7 +52,7 @@ "source": [ "qx = np.zeros_like(x)\n", "for i in range(len(x)):\n", - " qx[i], qy = ml.disvec(x[i], 1e-6)\n", + " qx[i : i + 1], qy = ml.disvec(x[i], 1e-6)\n", "plt.plot(x, qx)\n", "qxb = N * np.pi * R**2 / (2 * np.pi * R)\n", "plt.axhline(qxb, color=\"r\", ls=\"--\")\n", diff --git a/docs/steady/02examples/connected_wells.ipynb b/docs/steady/02examples/connected_wells.ipynb index 346b2a3..1f152b0 100755 --- a/docs/steady/02examples/connected_wells.ipynb +++ b/docs/steady/02examples/connected_wells.ipynb @@ -28,9 +28,19 @@ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", + "import timflow\n", "import timflow.steady as tfs" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "timflow.show_versions()" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -568,8 +578,7 @@ " levels=levels,\n", " decimals=2,\n", " layers=[ilay],\n", - " newfig=False,\n", - " linestyles=\"dashed\",\n", + " linestyles=\"solid\",\n", " color=\"C1\",\n", ")\n", "plt.plot(ws.xcp, ws.ycp, \"kx\");" @@ -718,7 +727,7 @@ " decimals=2,\n", " layers=[ilay],\n", " newfig=False,\n", - " linestyles=\"dashed\",\n", + " linestyles=\"solid\",\n", " color=\"C1\",\n", ")\n", "plt.plot(ws.xcp, ws.ycp, \"kx\");" diff --git a/docs/transient/00userguide/tutorials/tutorial0_start_a_model.ipynb b/docs/transient/00userguide/tutorials/tutorial0_start_a_model.ipynb index 82d7689..376bef9 100755 --- a/docs/transient/00userguide/tutorials/tutorial0_start_a_model.ipynb +++ b/docs/transient/00userguide/tutorials/tutorial0_start_a_model.ipynb @@ -35,7 +35,6 @@ "ls = tft.RiverString(ml, xy=[(-25, -30), (-15, 10), (10, 20)], tsandh=\"fixed\", layers=0)\n", "w = tft.Well(ml, xw=0, yw=0, rw=0.2, tsandQ=[(0, 1000)], layers=1)\n", "\n", - "ml.initialize()\n", "ml.solve(silent=True)" ] }, @@ -130,6 +129,13 @@ "source": [ "ax = ml.plots.head_along_line(-30, 10, 0, 0, t=10, layers=[0, 1, 2])" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/docs/transient/02examples/line_sink_well_sol.ipynb b/docs/transient/02examples/line_sink_well_sol.ipynb index 2fdb650..458a924 100755 --- a/docs/transient/02examples/line_sink_well_sol.ipynb +++ b/docs/transient/02examples/line_sink_well_sol.ipynb @@ -105,7 +105,7 @@ "metadata": {}, "outputs": [], "source": [ - "plt.figure(figsize=(16, 4))\n", + "plt.figure(figsize=(12, 4))\n", "for ilay in [0, 1, 2]:\n", " ax = plt.subplot(1, 3, ilay + 1)\n", " ml.plots.contour(\n", @@ -118,7 +118,8 @@ " labels=True,\n", " decimals=2,\n", " ax=ax,\n", - " )" + " )\n", + " ax.set_aspect(\"equal\", adjustable=\"box\")" ] }, { diff --git a/docs/transient/02examples/meandering_river.ipynb b/docs/transient/02examples/meandering_river.ipynb index 5c9f585..c5d158d 100644 --- a/docs/transient/02examples/meandering_river.ipynb +++ b/docs/transient/02examples/meandering_river.ipynb @@ -77,8 +77,22 @@ } ], "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, "language_info": { - "name": "python" + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.5" } }, "nbformat": 4, diff --git a/docs/transient/05benchmarks/well_benchmarks.ipynb b/docs/transient/05benchmarks/well_benchmarks.ipynb index 5e11a36..b5de1d6 100644 --- a/docs/transient/05benchmarks/well_benchmarks.ipynb +++ b/docs/transient/05benchmarks/well_benchmarks.ipynb @@ -181,7 +181,7 @@ "outputs": [], "source": [ "def volume(r, t=1):\n", - " return -2 * np.pi * r * ml.head(r, 0, t) * ml.aq.Scoefaq[0]\n", + " return (-2 * np.pi * r * ml.head(r, 0, t) * ml.aq.Scoefaq[0]).squeeze()\n", "\n", "\n", "quad(volume, 1e-5, np.inf)" diff --git a/pyproject.toml b/pyproject.toml index 813ba44..4ed63e6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -97,7 +97,7 @@ addopts = [ "--cov-report=term", "--cov-report=lcov:./coverage/lcov.info", "--ignore=docs/steady/02examples/vertical_anisotropy.ipynb", - "--ignore=docs/steady/04benchmarks/benchmarking_besselaes.ipynb", + "--ignore=docs/steady/04benchmarks/besselnumba_timing.ipynb", ] markers = [ "notebooks: run notebooks", diff --git a/timflow/steady/constant.py b/timflow/steady/constant.py index 4ad6db0..89d662b 100644 --- a/timflow/steady/constant.py +++ b/timflow/steady/constant.py @@ -170,7 +170,7 @@ def equation(self): # else: # mat[0, ieq:ieq+e. nunknowns] += -1 else: - rhs[0] -= e.potentiallayers( + rhs[0:1] -= e.potentiallayers( self.xc[icp], self.yc[icp], self.layers ).sum(0) return mat, rhs diff --git a/timflow/steady/linedoublet.py b/timflow/steady/linedoublet.py index 8d3ac31..36da118 100644 --- a/timflow/steady/linedoublet.py +++ b/timflow/steady/linedoublet.py @@ -408,8 +408,7 @@ def potinf(self, x, y, aq=None): rv = np.zeros((self.Nld, self.ldlist[0].nparam, aq.naq)) for i in range(self.Nld): rv[i] = self.ldlist[i].potinf(x, y, aq) - rv.shape = (self.nparam, aq.naq) - return rv + return rv.reshape((self.nparam, aq.naq)) def disvecinf(self, x, y, aq=None): if aq is None: @@ -417,8 +416,7 @@ def disvecinf(self, x, y, aq=None): rv = np.zeros((2, self.Nld, self.ldlist[0].nparam, aq.naq)) for i in range(self.Nld): rv[:, i] = self.ldlist[i].disvecinf(x, y, aq) - rv.shape = (2, self.nparam, aq.naq) - return rv + return rv.reshape((2, self.nparam, aq.naq)) def plot(self, ax=None, layer=None): if ax is None: diff --git a/timflow/steady/linesink.py b/timflow/steady/linesink.py index fea0539..9bfbf9a 100644 --- a/timflow/steady/linesink.py +++ b/timflow/steady/linesink.py @@ -597,7 +597,7 @@ def initialize(self): self.resfac = ( np.tile(self.res / self.whfac, self.ncp) * self.strengthinf ) # this is resfach ! - self.resfac.shape = (self.ncp, self.nlayers, self.nunknowns) + self.resfac = self.resfac.reshape((self.ncp, self.nlayers, self.nunknowns)) if len(self.hls) == 1: self.hc = self.hls * np.ones(self.nparam) elif len(self.hls) == self.ncp: # head specified at control points @@ -707,7 +707,7 @@ def equation(self): mat[0, ieq : ieq + e.nunknowns] = self.dischargeinf() break ieq += e.nunknowns - rhs[0] = self.Qls + rhs[0:1] = self.Qls return mat, rhs def setparams(self, sol): @@ -740,14 +740,16 @@ def __init__( self.xy = np.atleast_2d(xy).astype("d") if closed: self.xy = np.vstack((self.xy, self.xy[0])) - self.x, self.y = self.xy[:, 0], self.xy[:, 1] self.order = order # same for all segments in string self.dely = dely # same for all segments in string self.lslist = [] if self.xy.shape[1] == 2: self.nls = len(self.xy) - 1 + self.x, self.y = self.xy[:, 0], self.xy[:, 1] elif self.xy.shape[1] == 4: self.nls = len(self.xy) + self.x = None + self.y = None if self.layers.ndim == 1: if len(self.layers) == self.nls: self.layers = self.layers[:, np.newaxis] @@ -801,8 +803,7 @@ def potinf(self, x, y, aq=None): if aq in self.aq: for i, ls in enumerate(self.lslist): rv[i] = ls.potinf(x, y, aq) - rv.shape = (self.nparam, aq.naq) - return rv + return rv.reshape((self.nparam, aq.naq)) def disvecinf(self, x, y, aq=None): if aq is None: @@ -811,8 +812,7 @@ def disvecinf(self, x, y, aq=None): if aq in self.aq: for i, ls in enumerate(self.lslist): rv[:, i] = ls.disvecinf(x, y, aq) - rv.shape = (2, self.nparam, aq.naq) - return rv + return rv.reshape((2, self.nparam, aq.naq)) def dischargeinf(self): rv = np.zeros((self.nls, self.lslist[0].nparam)) @@ -829,7 +829,7 @@ def discharge_per_linesink(self): array of shape (nlay, nlinesinks) """ Qls = self.parameters[:, 0] * self.dischargeinf() - Qls.shape = (self.nls, self.nlayers, self.order + 1) + Qls = Qls.reshape((self.nls, self.nlayers, self.order + 1)) Qls = Qls.sum(axis=2) rv = np.zeros((self.model.aq.naq, self.nls)) for i, q in enumerate(Qls): @@ -840,7 +840,7 @@ def discharge(self): """Discharge of the element in each layer.""" rv = np.zeros(self.aq[0].naq) Qls = self.parameters[:, 0] * self.dischargeinf() - Qls.shape = (self.nls, self.nlayers, self.order + 1) + Qls = Qls.reshape((self.nls, self.nlayers, self.order + 1)) Qls = np.sum(Qls, 2) for i, q in enumerate(Qls): rv[self.layers[i]] += q @@ -947,7 +947,6 @@ def __init__( self.res = res self.wh = wh self.model.add_element(self) - # TO DO: TEST FOR DIFFERENT AQUIFERS AND LAYERS def initialize(self): if len(self.hls) == 1: # one value @@ -978,6 +977,8 @@ def initialize(self): y1, y2 = self.y[i : i + 2] elif self.xy.shape[1] == 4: x1, y1, x2, y2 = self.xy[i] + else: + raise ValueError("xy should have shape (n, 2) or (n, 4)") self.lslist.append( River( self.model, @@ -1018,7 +1019,7 @@ def equation(self): # include resistance by computing position of coefficients in matrix # and subtracting resistance terms iself = self.model.elementlist.index(self) - jcol = np.sum(e.nunknowns for e in self.model.elementlist[:iself]) + jcol = np.sum([e.nunknowns for e in self.model.elementlist[:iself]]).astype(int) irow = 0 for ls in self.lslist: for icp in range(ls.ncp): @@ -1114,7 +1115,7 @@ def equation(self): mat[0, ieq : ieq + self.nunknowns] = self.dischargeinf() break ieq += e.nunknowns - rhs[0] = self.Qls + rhs[0:1] = self.Qls return mat, rhs diff --git a/timflow/steady/well.py b/timflow/steady/well.py index 1e3e622..998f2b8 100644 --- a/timflow/steady/well.py +++ b/timflow/steady/well.py @@ -82,7 +82,7 @@ def initialize(self): self.parameters[:, 0] = self.Qw self.resfac = self.res / (2 * np.pi * self.rw * self.aq.Haq[self.layers]) self.resfac = self.resfac * np.identity(self.nlayers) - self.resfac.shape = ( + self.resfac = self.resfac.reshape( self.ncp, self.nlayers, self.nlayers, # changed to nlayers from nunknowns @@ -451,7 +451,7 @@ def equation(self): mat[0, ieq : ieq + e.nunknowns] = 1.0 break ieq += e.nunknowns - rhs[0] = self.Qw + rhs[0:1] = self.Qw return mat, rhs def setparams(self, sol): @@ -618,7 +618,7 @@ def equation(self): rhs[i] -= rhs[0] # first equation is head at control point equals hcp mat[0] = 0.0 - rhs[0] = self.hcp + rhs[0:1] = self.hcp aq = self.model.aq.find_aquifer_data(self.xcp, self.ycp) ieq = 0 for e in self.model.elementlist: @@ -628,7 +628,9 @@ def equation(self): ) ieq += e.nunknowns else: - rhs[0] -= e.potentiallayers(self.xcp, self.ycp, self.lcp) / aq.T[self.lcp] + rhs[0:1] -= ( + e.potentiallayers(self.xcp, self.ycp, self.lcp) / aq.T[self.lcp] + ) return mat, rhs def setparams(self, sol): @@ -727,7 +729,7 @@ def potinf(self, x, y, aq=None): if r < self.rw: r = self.rw # If at well, set to at radius if aq.ilap: - pot[0] = np.log(r / self.rw) / (2 * np.pi) + pot[0:1] = np.log(r / self.rw) / (2 * np.pi) pot[1:] = -k0(r / aq.lab[1:]) / (2 * np.pi) / k0(self.rw / aq.lab[1:]) else: pot[:] = -k0(r / aq.lab) / (2 * np.pi) / k0(self.rw / aq.lab) @@ -751,8 +753,8 @@ def disvecinf(self, x, y, aq=None): xminxw = self.rw yminyw = 0.0 if aq.ilap: - qx[0] = -1 / (2 * np.pi) * xminxw / rsq - qy[0] = -1 / (2 * np.pi) * yminyw / rsq + qx[0:1] = -1 / (2 * np.pi) * xminxw / rsq + qy[0:1] = -1 / (2 * np.pi) * yminyw / rsq kone = k1(r / aq.lab[1:]) / k0(self.rw / aq.lab[1:]) qx[1:] = -kone * xminxw / (r * aq.lab[1:]) / (2 * np.pi) qy[1:] = -kone * yminyw / (r * aq.lab[1:]) / (2 * np.pi) @@ -846,7 +848,7 @@ def potinf(self, x, y, aq=None): for w in self.wlist: rv[j : j + w.nparam] = w.potinf(x, y, aq) j += w.nparam - # rv.shape = (self.nparam, aq.naq) + # rv = rv.reshape((self.nparam, aq.naq)) return rv def disvecinf(self, x, y, aq=None): @@ -907,7 +909,7 @@ def equation(self): # include resistance by finding position of element in coefficient matrix # and subtracting resfac (resistance factor). iself = self.model.elementlist.index(self) - jcol = np.sum(e.nunknowns for e in self.model.elementlist[:iself]) + jcol = np.sum([e.nunknowns for e in self.model.elementlist[:iself]]) irow = 0 for w in self.wlist: mat[irow : irow + w.nlayers, jcol : jcol + w.nunknowns] -= w.resfac[ @@ -1004,7 +1006,7 @@ def equation(self): mat[0, ieq : ieq + self.nunknowns] = 1.0 break ieq += e.nunknowns - rhs[0] = self.Qw + rhs[0:1] = self.Qw return mat, rhs def setparams(self, sol): @@ -1166,7 +1168,7 @@ def equation(self): # first equation is head at control point equals hcp mat[0] = 0 - rhs[0] = self.hcp + rhs[0:1] = self.hcp aq = self.model.aq.find_aquifer_data(self.xcp, self.ycp) ieq = 0 for e in self.model.elementlist: @@ -1184,7 +1186,9 @@ def equation(self): ) ieq += e.nunknowns else: - rhs[0] -= e.potentiallayers(self.xcp, self.ycp, self.lcp) / aq.T[self.lcp] + rhs[0:1] -= ( + e.potentiallayers(self.xcp, self.ycp, self.lcp) / aq.T[self.lcp] + ) return mat, rhs def setparams(self, sol): diff --git a/timflow/transient/aquifer.py b/timflow/transient/aquifer.py index 4013083..a7e5016 100644 --- a/timflow/transient/aquifer.py +++ b/timflow/transient/aquifer.py @@ -135,10 +135,12 @@ def initialize(self): self.coef[:, :, i] = np.linalg.solve(v, b).T self.lab = 1.0 / np.sqrt(self.eigval) self.lab2 = self.lab.copy() - self.lab2.shape = (self.naq, self.model.nint, self.model.npint) + self.lab2 = self.lab2.reshape((self.naq, self.model.nint, self.model.npint)) self.lababs = np.abs(self.lab2[:, :, 0]) # used to check distances self.eigvec2 = self.eigvec.copy() - self.eigvec2.shape = (self.naq, self.naq, self.model.nint, self.model.npint) + self.eigvec2 = self.eigvec2.reshape( + (self.naq, self.naq, self.model.nint, self.model.npint) + ) def compute_lab_eigvec(self, p, returnA=False, B=None): sqrtpSc = np.sqrt(p * self.Scoefll * self.c) diff --git a/timflow/transient/circareasink.py b/timflow/transient/circareasink.py index 1ac364a..83797c2 100644 --- a/timflow/transient/circareasink.py +++ b/timflow/transient/circareasink.py @@ -66,7 +66,7 @@ def initialize(self): self.setflowcoef() # Since recharge is in layer 0, and RHS is -N self.an = self.aq.coef[0, :] * self.flowcoef - self.an.shape = (self.aq.naq, self.model.nint, self.model.npint) + self.an = self.an.reshape((self.aq.naq, self.model.nint, self.model.npint)) self.termin = self.aq.lab2 * self.R * self.an self.termin2 = self.aq.lab2**2 * self.an self.terminq = self.R * self.an @@ -105,8 +105,7 @@ def potinf(self, x, y, aq=None): for j in range(self.model.nint): if (r - self.R) / abs(self.aq.lab2[i, j, 0]) < self.rzero: rv[0, i, j, :] = self.termout[i, j, :] * self.I1RK0r(r, i, j) - rv.shape = (self.nparam, aq.naq, self.model.npval) - return rv + return rv.reshape((self.nparam, aq.naq, self.model.npval)) def disvecinf(self, x, y, aq=None): """Can be called with only one x,y value.""" @@ -129,7 +128,7 @@ def disvecinf(self, x, y, aq=None): for j in range(self.model.nint): if (r - self.R) / abs(self.aq.lab2[i, j, 0]) < self.rzero: qr[0, i, j] = self.termoutq[i, j, :] * self.I1RK1r(r, i, j) - qr.shape = (self.nparam, aq.naq, self.model.npval) + qr = qr.reshape((self.nparam, aq.naq, self.model.npval)) qx[:] = qr * (x - self.xc) / r qy[:] = qr * (y - self.yc) / r return qx, qy diff --git a/timflow/transient/circinhom.py b/timflow/transient/circinhom.py index 10ab21d..9a2bf4b 100644 --- a/timflow/transient/circinhom.py +++ b/timflow/transient/circinhom.py @@ -274,8 +274,7 @@ def potinf(self, x, y, aq=None): rv[self.aqin.Naq + i, i, j, :] = self.approx.kvratio( r, self.R, self.aqout.lab2[i, j, :] ) - rv.shape = (self.Nparam, aq.Naq, self.model.Np) - return rv + return rv.reshape((self.Nparam, aq.Naq, self.model.Np)) def disinf(self, x, y, aq=None): """Can be called with only one x,y value.""" @@ -304,7 +303,7 @@ def disinf(self, x, y, aq=None): -self.approx.ivratiop(r, self.R, self.aqin.lab2[i, j, :]) / self.aqin.lab2[i, j, :] ) - qr.shape = (self.nparam, aq.naq, self.model.np) + qr = qr.reshape((self.nparam, aq.naq, self.model.np)) qx[:] = qr * (x - self.x0) / r qy[:] = qr * (y - self.y0) / r if aq == self.aqout: @@ -326,7 +325,7 @@ def disinf(self, x, y, aq=None): self.approx.kvratiop(r, self.R, self.aqout.lab2[i, j, :]) / self.aqout.lab2[i, j, :] ) - qr.shape = (self.Nparam, aq.Naq, self.model.Np) + qr = qr.reshape((self.Nparam, aq.Naq, self.model.Np)) qx[:] = qr * (x - self.x0) / r qy[:] = qr * (y - self.y0) / r return qx, qy @@ -427,7 +426,7 @@ def layout(self): # for n in range(1,self.order+1): # rv[aq.Naq+i,2*n-1,i,j,:] = pot[n] * np.cos(n*alpha) # rv[aq.Naq+i,2*n ,i,j,:] = pot[n] * np.sin(n*alpha) -# rv.shape = (self.Nparam,aq.Naq,self.model.Np) +# rv = rv.reshape((self.Nparam, aq.Naq, self.model.Np)) # return rv # def disinf(self,x,y,aq=None): # '''Can be called with only one x,y value''' @@ -462,8 +461,8 @@ def layout(self): # qr[i,2*n ,i,j,:] = -potp[n] / 2 / self.aqin.lab2[i,j,:] * np.sin(n*alpha) # qt[i,2*n-1,i,j,:] = pot[n] * np.sin(n*alpha) * n / r # qt[i,2*n ,i,j,:] = -pot[n] * np.cos(n*alpha) * n / r -# qr.shape = (self.Nparam/2,aq.Naq,self.model.Np) -# qt.shape = (self.Nparam/2,aq.Naq,self.model.Np) +# qr = qr.reshape((self.Nparam / 2, aq.Naq, self.model.Np)) +# qt = qt.reshape((self.Nparam / 2, aq.Naq, self.model.Np)) # qx[:self.Nparam/2,:,:] = qr * np.cos(alpha) - qt * np.sin(alpha); # qy[:self.Nparam/2,:,:] = qr * np.sin(alpha) + qt * np.cos(alpha); # if aq == self.aqout: @@ -494,8 +493,8 @@ def layout(self): # qr[i,2*n ,i,j,:] = -potp[n] / self.aqout.lab2[i,j,:] * np.sin(n*alpha) # qt[i,2*n-1,i,j,:] = pot[n] * np.sin(n*alpha) * n / r # qt[i,2*n ,i,j,:] = -pot[n] * np.cos(n*alpha) * n / r -# qr.shape = (self.Nparam/2,aq.Naq,self.model.Np) -# qt.shape = (self.Nparam/2,aq.Naq,self.model.Np) +# qr = qr.reshape((self.Nparam / 2, aq.Naq, self.model.Np)) +# qt = qt.reshape((self.Nparam / 2, aq.Naq, self.model.Np)) # qx[self.Nparam/2:,:,:] = qr * np.cos(alpha) - qt * np.sin(alpha); # qy[self.Nparam/2:,:,:] = qr * np.sin(alpha) + qt * np.cos(alpha); # return qx,qy diff --git a/timflow/transient/linedoublet.py b/timflow/transient/linedoublet.py index 4bbf0e1..29c2234 100644 --- a/timflow/transient/linedoublet.py +++ b/timflow/transient/linedoublet.py @@ -143,8 +143,7 @@ def potinf(self, x, y, aq=None): ) for k in range(self.nlayers): rv[k :: self.nlayers, i, j, :] = self.term2[k, i, j, :] * pot - rv.shape = (self.nparam, aq.naq, self.model.npval) - return rv + return rv.reshape((self.nparam, aq.naq, self.model.npval)) def disvecinf(self, x, y, aq=None): """Can be called with only one x,y value.""" @@ -182,10 +181,10 @@ def disvecinf(self, x, y, aq=None): rvy[k :: self.nlayers, i, j, :] = ( self.term2[k, i, j, :] * qxqy[self.order + 1 :, :] ) - - rvx.shape = (self.nparam, aq.naq, self.model.npval) - rvy.shape = (self.nparam, aq.naq, self.model.npval) - return rvx, rvy + return ( + rvx.reshape((self.nparam, aq.naq, self.model.npval)), + rvy.reshape((self.nparam, aq.naq, self.model.npval)), + ) def plot(self, ax=None, layer=None): if ax is None: diff --git a/timflow/transient/linesink.py b/timflow/transient/linesink.py index 06f8c24..0109a6a 100644 --- a/timflow/transient/linesink.py +++ b/timflow/transient/linesink.py @@ -138,8 +138,7 @@ def potinf(self, x, y, aq=None): ) # Divide by L as the parameter is total discharge rv[:, i, j, :] = self.term2[:, i, j, :] * pot / self.L - rv.shape = (self.nparam, aq.naq, self.model.npval) - return rv + return rv.reshape((self.nparam, aq.naq, self.model.npval)) def disvecinf(self, x, y, aq=None): """Can be called with only one x,y value.""" @@ -172,9 +171,10 @@ def disvecinf(self, x, y, aq=None): ) rvx[:, i, j, :] = self.term2[:, i, j, :] * qxqy[0] rvy[:, i, j, :] = self.term2[:, i, j, :] * qxqy[1] - rvx.shape = (self.nparam, aq.naq, self.model.npval) - rvy.shape = (self.nparam, aq.naq, self.model.npval) - return rvx, rvy + return ( + rvx.reshape((self.nparam, aq.naq, self.model.npval)), + rvy.reshape((self.nparam, aq.naq, self.model.npval)), + ) def headinside(self, t): """The head inside the line-sink. @@ -419,8 +419,8 @@ def initialize(self): self.dischargeinflayers[i * self.nlayers : (i + 1) * self.nlayers, :] = ( self.lslist[i].dischargeinflayers ) - self.xc[i] = self.lslist[i].xc - self.yc[i] = self.lslist[i].yc + self.xc[i : i + 1] = self.lslist[i].xc + self.yc[i : i + 1] = self.lslist[i].yc def potinf(self, x, y, aq=None): """Returns array (nunknowns, Nperiods).""" @@ -976,9 +976,7 @@ def potinf(self, x, y, aq=None): ) for k in range(self.nlayers): rv[k :: self.nlayers, i, j, :] = self.term2[k, i, j, :] * pot - - rv.shape = (self.nparam, aq.naq, self.model.npval) - return rv + return rv.reshape((self.nparam, aq.naq, self.model.npval)) def disvecinf(self, x, y, aq=None): """Can be called with only one x,y value.""" @@ -1016,9 +1014,10 @@ def disvecinf(self, x, y, aq=None): rvy[k :: self.nlayers, i, j, :] = ( self.term2[k, i, j, :] * qxqy[self.order + 1 :, :] ) - rvx.shape = (self.nparam, aq.naq, self.model.npval) - rvy.shape = (self.nparam, aq.naq, self.model.npval) - return rvx, rvy + return ( + rvx.reshape((self.nparam, aq.naq, self.model.npval)), + rvy.reshape((self.nparam, aq.naq, self.model.npval)), + ) def headinside(self, t): """The head inside the line-sink. diff --git a/timflow/transient/model.py b/timflow/transient/model.py index b56d2ae..83801a3 100644 --- a/timflow/transient/model.py +++ b/timflow/transient/model.py @@ -416,9 +416,9 @@ def velocomp(self, x, y, z, t, aq=None, layer_ltype=None): )[:, 0] else: # layer = 0, so top layer if aq.naq == 1: # only one layer - h[1] = self.head(x, y, t, layers=[layer], aq=aq, neglect_steady=True)[ - :, 0 - ] + h[1:2] = self.head( + x, y, t, layers=[layer], aq=aq, neglect_steady=True + )[:, 0] else: h[1:] = self.head( x, y, t, layers=[layer, layer + 1], aq=aq, neglect_steady=True diff --git a/timflow/transient/well.py b/timflow/transient/well.py index e819cb5..5d5a8ee 100644 --- a/timflow/transient/well.py +++ b/timflow/transient/well.py @@ -109,8 +109,7 @@ def potinf(self, x, y, aq=None): # quicker? # bessel.k0besselv( r / self.aq.lab2[i,j,:], pot ) rv[:, i, j, :] = self.term2[:, i, j, :] * pot - rv.shape = (self.nparam, aq.naq, self.model.npval) - return rv + return rv.reshape((self.nparam, aq.naq, self.model.npval)) def potinfone(self, x, y, jtime, aq=None): """Can be called with only one x,y value for time interval jtime.""" @@ -126,7 +125,7 @@ def potinfone(self, x, y, jtime, aq=None): if r / abs(self.aq.lab2[i, jtime, 0]) < self.rzero: pot[:] = kv(0, r / self.aq.lab2[i, jtime, :]) rv[:, i, :] = self.term2[:, i, jtime, :] * pot - # rv.shape = (self.nparam, aq.naq, self.model.npval) + # rv = rv.reshape((self.nparam, aq.naq, self.model.npval)) return rv def disvecinf(self, x, y, aq=None): @@ -151,7 +150,7 @@ def disvecinf(self, x, y, aq=None): * kv(1, r / self.aq.lab2[i, j, :]) / self.aq.lab2[i, j, :] ) - qr.shape = (self.nparam, aq.naq, self.model.npval) + qr = qr.reshape((self.nparam, aq.naq, self.model.npval)) qx[:] = qr * (x - self.xw) / r qy[:] = qr * (y - self.yw) / r return qx, qy