-
Couldn't load subscription status.
- Fork 5
Neuman sloution Moench et al. 2001 #6
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
base: main
Are you sure you want to change the base?
Conversation
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It's a starting point. 😉
anaflow/flow/Neuman.py
Outdated
| __all__ = [] | ||
|
|
||
|
|
||
| def neuman( |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
rename to neuman_unconfined_laplace
anaflow/flow/Neuman.py
Outdated
| return res | ||
|
|
||
|
|
||
| def neuman_from_laplace( |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
rename to neuman_unconfined
anaflow/flow/Neuman.py
Outdated
| Parameters | ||
| ---------- | ||
| s : :class:`numpy.ndarray` | ||
| Array with all time-points where the function should be evaluated in the Laplace space. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
"Laplace-space-points" not "time-points"
anaflow/flow/Neuman.py
Outdated
| kz : :class:`float`, optional | ||
| Hydraulic conductivity in the vertical direction | ||
| kr : :class:`float`, optional | ||
| Hydraulic conductivity in the horizontal direction |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
actually we don't need these two values but rather a single anis (in line with other solutions) value for the anisotropy ratio kz/kr. This is named kd in the Moench paper.
anaflow/flow/Neuman.py
Outdated
| l : :class:`float`, optional | ||
| Vertical distance from initial water table to bottom | ||
| of pumped well screen |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
would rename to well_depth to be in line with welltestpy.
anaflow/flow/Neuman.py
Outdated
| d : :class:`float`, optional | ||
| Vertical distance from initial water table to top of | ||
| pumped well screen |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Would use screen_size instead and calculate d = well_depth - screen_size to be in line with GeoStat-Framework/welltestpy#18
anaflow/flow/Neuman.py
Outdated
| rw : :class:`float`, optional | ||
| Outside radius of the pumped well screen |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I would use the solution from Moench 1993, where infinitely small wells are assumed. I think we will run into numerical problems, that we don't have the time to solve. So I would vote for leaving rw out for now.
anaflow/flow/Neuman.py
Outdated
| Hydraulic conductivity in the vertical direction | ||
| kr : :class:`float`, optional | ||
| Hydraulic conductivity in the horizontal direction | ||
| Ss : :class:`float`, optional |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
would rename to specific_storage to be more explicit.
anaflow/flow/Neuman.py
Outdated
| time = [10, 600, 36000] # 10s, 10min, 10h | ||
| rad = np.geomspace(0.1, 10) | ||
| # default parameters | ||
| T = 1e-4 | ||
| S = 1e-4 | ||
| Q = -1e-4 | ||
|
|
||
| neu = neuman_from_laplace(time, rad, storage=S, transmissivity=T) | ||
| for i in range(len(time)): | ||
| plt.plot(rad, neu[i, :], color="C0", linestyle=":", linewidth=4) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Don't forget to remove this after testing.
|
Here a proposal for an even better root finder for from scipy.optimize import root_scalar
def get_f_df_df2(value=0):
"""Get target function and its first two derivatives."""
if value < 0:
raise ValueError("Neuman: epsilon for root finding needs to be >0.")
def _f_df_df2(x):
return (
np.subtract(np.multiply(x, np.tan(x)), value),
np.tan(x) + np.divide(x, np.cos(x) ** 2),
2 * (np.multiply(x, np.tan(x)) + 1.0) * np.cos(x) ** -2,
)
return _f_df_df2
def nth_root(n, v):
"""Get n-th root of x*tan(x) = v."""
x0 = np.sqrt(v) if (v < 1 and n < 1) else np.arctan(v) + n * np.pi
f = get_f_df_df2(v)
sol = root_scalar(f, x0=x0, fprime=True, fprime2=True)
if not sol.converged:
raise ValueError(f"Neuman: couldn't find {n}-th root for eps={v}")
return sol.root |
|
@JarnoHerr do you still have interest and/or time to finish this, or should I take over? Greetings! |
|
@MuellerSeb the interest is still there, but time is a bit limited at this time and I really don't want to halt the progress of the package, so you can take over. I look forward to the progress that will be made. Best Regards! |
No description provided.