Skip to content
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

Neuman sloution Moench et al. 2001 #6

Draft
wants to merge 9 commits into
base: main
Choose a base branch
from

Conversation

JarnoHerr
Copy link

No description provided.

@MuellerSeb MuellerSeb self-requested a review June 16, 2021 13:34
@MuellerSeb MuellerSeb added the enhancement New feature or request label Jun 16, 2021
@MuellerSeb MuellerSeb added this to the v1.1 milestone Jun 16, 2021
Copy link
Member

@MuellerSeb MuellerSeb left a 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. 😉

__all__ = []


def neuman(
Copy link
Member

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

return res


def neuman_from_laplace(
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

rename to neuman_unconfined

Parameters
----------
s : :class:`numpy.ndarray`
Array with all time-points where the function should be evaluated in the Laplace space.
Copy link
Member

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"

Comment on lines 67 to 70
kz : :class:`float`, optional
Hydraulic conductivity in the vertical direction
kr : :class:`float`, optional
Hydraulic conductivity in the horizontal direction
Copy link
Member

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.

Comment on lines 62 to 64
l : :class:`float`, optional
Vertical distance from initial water table to bottom
of pumped well screen
Copy link
Member

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.

Comment on lines 59 to 61
d : :class:`float`, optional
Vertical distance from initial water table to top of
pumped well screen
Copy link
Member

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

Comment on lines 65 to 66
rw : :class:`float`, optional
Outside radius of the pumped well screen
Copy link
Member

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.

Hydraulic conductivity in the vertical direction
kr : :class:`float`, optional
Hydraulic conductivity in the horizontal direction
Ss : :class:`float`, optional
Copy link
Member

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.

Comment on lines 153 to 162
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)
Copy link
Member

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.

@MuellerSeb
Copy link
Member

Here a proposal for an even better root finder for x tan(x) that uses the second derivative:

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

@MuellerSeb
Copy link
Member

@JarnoHerr do you still have interest and/or time to finish this, or should I take over? Greetings!

@JarnoHerr
Copy link
Author

@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!

@MuellerSeb MuellerSeb changed the base branch from develop to main April 16, 2023 10:11
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants