diff --git a/tests/test_connector.py b/tests/test_connector.py index cbadb4b..ac94dac 100644 --- a/tests/test_connector.py +++ b/tests/test_connector.py @@ -1,7 +1,8 @@ import numpy as np -from welleng.connector import Connector +from welleng.connector import Connector, drop_off, extend_to_tvd from welleng.survey import Survey, from_connections +from welleng.node import Node def test_md_hold(): @@ -111,3 +112,65 @@ def test_radius_critical_with_min_curve(): azi2=0, ) assert c.radius_critical < c.radius_design + +def test_drop_off(tol=1e-4): + node = Node( + pos=[0, 0, 3000], inc=30, azi=135, md=4000 + ) + nodes = drop_off( + target_inc=0, dls=3, node=node + ) + assert len(nodes) == 1, "Unexpected tangent section." + assert abs(nodes[0].inc_deg - 0) < tol + + nodes = drop_off( + target_inc=0, dls=3, node=node, delta_md=1000 + ) + assert len(nodes) == 2, "Unexpected number of nodes." + assert nodes[-1].md - 1000 - node.md < tol, "Unexpected delta md." + assert nodes[0].inc_deg - nodes[1].inc_deg < tol, "Unexpected tangent inc." + +def test_extend_to_tvd(tol=1e-4): + node = Node( + pos=[0, 0, 3000], inc=30, azi=135, md=4000 + ) + nodes = drop_off( + target_inc=0, dls=3, node=node + ) + connectors = extend_to_tvd( + target_tvd=3500, node=node, target_inc=0, dls=3 + ) + assert len(connectors) == 2, "Unexpected number of connectors." + assert connectors[0].node_end.md - nodes[0].md < tol, "Unexpected md." + assert np.allclose( + connectors[0].node_end.pos_nev, + nodes[0].pos_nev + ), "Unexpected pos_nev." + assert np.allclose( + connectors[0].node_end.vec_nev, + nodes[0].vec_nev, rtol=tol, atol=tol + ), "Unexpected vec_nev." + assert np.allclose( + connectors[0].node_end.vec_nev, + connectors[1].node_end.vec_nev, rtol=tol, atol=tol + ), "Unexpected tangent section." + assert connectors[1].node_end.pos_nev[2] - 3500 < tol, "Unexpected tvd." + + connectors = extend_to_tvd( + target_tvd=3500, node=node + ) + assert connectors[-1].node_end.pos_nev[2] - 3500 < tol, "Unexpected tvd." + assert np.allclose( + connectors[-1].node_start.vec_nev, + connectors[-1].node_end.vec_nev, rtol=tol, atol=tol + ), "Unexpected tangent section." + pass + + +def main(): + test_drop_off() + test_extend_to_tvd() + + +if __name__ == "__main__": + main() diff --git a/tests/test_utils.py b/tests/test_utils.py index 3c726fd..70f7021 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -42,7 +42,7 @@ def _generate_random_dms(n: int, ndigits: int = None) -> NDArray: ] ) - return np.stack((deg, min, sec, direction), axis=-1, dtype=object).reshape( + return np.stack((deg, min, sec, direction), axis=-1, dtype='object').reshape( (-1, 2, 4) ) @@ -232,3 +232,12 @@ def test_get_toolface(): pos2 = pos1 + np.array([radius, 0, radius]) toolface = get_toolface(pos1, vec1, pos2) assert np.isclose(toolface, 0) + + +def main(): + test_dms2decimal2dms() + pass + + +if __name__ == "__main__": + main() diff --git a/welleng/connector.py b/welleng/connector.py index 44c1138..e5ca3ed 100644 --- a/welleng/connector.py +++ b/welleng/connector.py @@ -1,7 +1,7 @@ from copy import copy, deepcopy import numpy as np -# from scipy.optimize import minimize +from scipy.optimize import minimize from scipy.spatial import distance try: @@ -11,8 +11,11 @@ NUMBA = False from .node import Node, get_node_params -from .utils import (NEV_to_HLA, _get_angles, dls_from_radius, get_angles, - get_nev, get_unit_vec, get_vec, get_xyz) +from .utils import ( + NEV_to_HLA, _get_angles, dls_from_radius, get_angles, + get_nev, get_unit_vec, get_vec, get_xyz, radius_from_dls, + get_arc +) class Connector: @@ -1632,6 +1635,237 @@ def _get_section( ) +def drop_off( + target_inc: float, dls: float, delta_md: float | None = None, + node: Node | None = None, tol: float = 1e-5 +) -> list: + """Determines the sections necessary to drop off (or build) towards a + target inclination, constrained by the section length and dls. + + Recommended to use the ``extend_to_tvd`` function if a specific tvd + target is required. + + Parameters + ---------- + target_inc: float + The target trajectory inclination in degrees. + dls: float + The design dls in deg/30m. + delta_md: float + The section length in meters. + node: ``Node`` + The ``Node`` of the current location. + + Returns + ------- + nodes: list + A list of ``Connectors`` that describe the required trajectory changes. + The list contains at least one value (describing the arc or curve) and + a second value if the target inclination was achieved, to maintain a + tangent section. + """ + def _drop_off( + x: tuple, + node, + return_data: bool = False + ) -> tuple: + dogleg, toolface = x + pos2, vec2, arc_length = get_arc( + dogleg, radius, + 0 if -np.pi / 2 <= toolface <= np.pi / 2 else np.pi, + node.pos_nev, node.vec_nev + ) + if return_data: + return (pos2, vec2, arc_length) + else: + inc, _ = np.degrees( + get_angles(vec2, nev=True)[0] + ) + return abs(inc - target_inc) + + node = Node() if node is None else node + node.md = 0 if node.md is None else node.md + radius = radius_from_dls(dls) + if isinstance(delta_md, np.ndarray): + delta_md = delta_md[0] + max_dogleg = ( + 2 * np.pi if delta_md is None + else delta_md / radius + ) + + args = (node,) + bounds = [ + [0, min(2 * np.pi, max_dogleg)], + [-np.pi, np.pi] + ] + x0 = [bounds[0][1] / 2, np.pi] + result = minimize( + _drop_off, x0, args=args, bounds=bounds, + method="SLSQP" + ) + pos2, vec2, arc_length = _drop_off( + result.x, *args, return_data=True + ) + tangent_length = ( + 0 if delta_md is None + else delta_md - arc_length + ) + node2 = Node( + pos=pos2, vec=vec2, md=node.md + arc_length + ) + if tangent_length > tol: + pos3 = pos2 + tangent_length * vec2 + node3 = Node( + pos=pos3, vec=vec2, md=node2.md + tangent_length + ) + return [node2, node3] + else: + return [node2] + + +def extend_to_tvd( + target_tvd: float, node: Node | None = None, + delta_md: float | None = None, + target_inc: float | None = None, dls: float | None = None +) -> list: + """Determines a list of ``Connector`` sections required to drop off (or + build) to achieve a target tvd. + + Parameters + ---------- + target_tvd: float + The target true vertical depth in meters. + node: ``Node`` | None (default=None) + The ``Node`` describing the current well position. If not provided, + defaults to surface and pointing down. + delta_md: float | None (default=None) + Will limit the section length to the value provided in meters. + target_inc: float | None (default=None) + The target inclination at the target tvd in degrees. If provided, + will get as close to this value as feasible with the provided dls + constraint. If the target inclination can be achieved, a tangent + section will be held until the target tvd is achieved. + dls: float | None (default=None) + The design dls in deg/30m. If none provided, defaults to 2.3 deg/30m. + + Returns + ------- + nodes: list + A list of ``Connector`` that describe the required trajectory changes. + The list contains at least one value (describing the arc or curve) and + a second value if the target inclination was achieved, to maintain a + tangent section. + + Examples + -------- + A well has landed at the top reservoir with an inclination of 30 degrees + and the Wells Engineer wants to drop this angle to 0 degrees while drilling + the reservoir section. + + >>> import welleng as we + >>> from pprint import pprint + >>> node = we.node.Node(pos=[0, 0, 3000], md=4000, inc=30, azi=135) + >>> connectors = we.connector.extend_to_tvd(target_tvd=3200, node=node, target_inc=0, dls=3) + >>> len(connectors) + ... 1 + + Since there's only one ``Connector`` in ``connectors`` then likely the + target inclination has not been met, else there would also be a tangent + section in this case. + + >>> pprint(connectors[-1].node_end.__dict__) + ... {'azi_deg': 134.99999999999994, + 'azi_rad': 2.356194490192344, + 'cov_nev': array([[[0., 0., 0.], + [0., 0., 0.], + [0., 0., 0.]]]), + 'inc_deg': 8.681065764308594, + 'inc_rad': 0.151513180169343, + 'md': 4213.189342356914, + 'pos_nev': [-49.63739779370637, 49.63739779370647, 3199.9999999980705], + 'pos_xyz': [49.63739779370647, -49.63739779370637, 3199.9999999980705], + 'unit': 'meters', + 'vec_nev': [-0.10672656069796796, 0.1067265606979681, 0.9885438192023488], + 'vec_xyz': [0.1067265606979681, -0.10672656069796796, 0.9885438192023488]} + + Indeed, the ``inc_deg`` attribute shows the inclination dropped to just + below 9 degrees and the z coordinate in ``pos_nev`` shows that the target + tvd was met. + """ + if node is None: + node = Node() + node.md = 0 if node.md is None else node.md # default value is None which complicates things + _delta_md = 1e-8 if delta_md is None else delta_md + connections = [] + if target_inc is None: + def _extend_tvd( + delta_md, pos, vec, target_tvd, return_data=False + ): + pos2 = np.array(pos) + delta_md * np.array(vec) + if return_data: + return pos2 + else: + return abs(pos2[2] - target_tvd) + + args = (node.pos_nev, node.vec_nev, target_tvd) + bounds = [[0, None]] + result = minimize( + _extend_tvd, + _delta_md, + args=args, + bounds=bounds, + method="Powell" + ) + pos2 = _extend_tvd(result.x, *args, return_data=True) + connections.append(Connector( + pos1=node.pos_nev, vec1=node.vec_nev, + md1=0 if node.md is None else node.md, + pos2=pos2, + dls_design=dls, force_min_curve=True + )) + else: + def _drop_off( + delta_md, target_inc, target_tvd, dls, node, + return_data=False + ): + nodes = drop_off( + target_inc, dls, delta_md, node + ) + if return_data: + return nodes + else: + return abs(nodes[-1].pos_nev[2] - target_tvd) + + dls = 2.5 if dls is None else dls + args = (target_inc, target_tvd, dls, node) + bounds = [ + [min(target_tvd - node.pos_nev[2], _delta_md), None] + ] + x0 = _delta_md + result = minimize( + _drop_off, + x0, + args=args, + bounds=bounds, + method='Powell' + ) + nodes = _drop_off( + result.x, *args, return_data=True + ) + connections.append(Connector( + node1=node, + pos2=nodes[0].pos_nev, vec2=nodes[0].vec_nev, + dls_design=dls, force_min_curve=True + )) + if len(nodes) > 1: + connections.append(Connector( + node1=connections[-1].node_end, + pos2=nodes[-1].pos_nev, + dls_design=dls, force_min_curve=True + )) + return connections + + def numbafy(functions): for f in functions: f = njit(f) diff --git a/welleng/survey.py b/welleng/survey.py index 9de4284..a048c13 100644 --- a/welleng/survey.py +++ b/welleng/survey.py @@ -2247,14 +2247,14 @@ def func(x0, survey, dls_cont, tolerance): return diff -def _remove_duplicates(md, inc, azi): +def _remove_duplicates(md, inc, azi, decimals=4): arr = np.array([md, inc, azi]).T upper = arr[:-1] lower = arr[1:] temp = np.vstack(( upper[0], - lower[lower[:, 0] != upper[:, 0]] + lower[lower[:, 0].round(decimals) != upper[:, 0].round(decimals)] )) return temp.T @@ -2266,26 +2266,32 @@ def from_connections( start_xyz=[0., 0., 0.], start_cov_nev=None, radius=10, deg=False, error_model=None, - depth_unit='meters', surface_unit='meters' + depth_unit='meters', surface_unit='meters', + decimals: int | None = None ): """ Constructs a well survey from a list of sections of control points. Parameters ---------- - section_data: list of dicts with section data - start_nev: (3) array of floats (default: [0,0,0]) - The starting position in NEV coordinates. - radius: float (default: 10) - The radius is passed to the `welleng.survey.Survey` object - and represents the radius of the wellbore. It is also used - when visualizing the results, so can be used to make the - wellbore *thicker* in the plot. + section_data: list of dicts with section data + start_nev: (3) array of floats (default: [0,0,0]) + The starting position in NEV coordinates. + radius: float (default: 10) + The radius is passed to the `welleng.survey.Survey` object + and represents the radius of the wellbore. It is also used + when visualizing the results, so can be used to make the + wellbore *thicker* in the plot. + decimals: int (default=6) + Round the md decimal when checking for duplicate surveys. Results ------- survey: `welleng.survey.Survey` object """ + decimals = 6 if decimals is None else decimals + assert isinstance(decimals, int), "decimals must be an int" + if type(section_data) is not list: section_data = [section_data] diff --git a/welleng/utils.py b/welleng/utils.py index 2cfc196..6455ef4 100644 --- a/welleng/utils.py +++ b/welleng/utils.py @@ -497,6 +497,26 @@ def errors_from_cov(cov, data=False): return np.array([nn, ne, nv, ee, ev, vv]).T +def _get_arc_pos_and_vec(dogleg, radius): + pos = np.array([ + np.cos(dogleg), + 0., + np.sin(dogleg) + ]) * radius + pos[0] = radius - pos[0] + + vec = np.array([ + np.sin(dogleg), + 0., + np.cos(dogleg) + ]) + return (pos, vec) + + +if NUMBA: + _get_arc_pos_and_vec = njit(_get_arc_pos_and_vec) + + class Arc: def __init__(self, dogleg, radius): """ @@ -520,18 +540,20 @@ def __init__(self, dogleg, radius): self.radius = radius self.delta_md = dogleg * radius - self.pos = np.array([ - np.cos(dogleg), - 0., - np.sin(dogleg) - ]) * radius - self.pos[0] = radius - self.pos[0] - - self.vec = np.array([ - np.sin(dogleg), - 0., - np.cos(dogleg) - ]) + self.pos, self.vec = _get_arc_pos_and_vec(dogleg, radius) + + # self.pos = np.array([ + # np.cos(dogleg), + # 0., + # np.sin(dogleg) + # ]) * radius + # self.pos[0] = radius - self.pos[0] + + # self.vec = np.array([ + # np.sin(dogleg), + # 0., + # np.cos(dogleg) + # ]) def transform(self, toolface, pos=None, vec=None, target=False): """ @@ -579,10 +601,11 @@ def transform(self, toolface, pos=None, vec=None, target=False): return (pos_new, vec_new) -def get_arc(dogleg, radius, toolface, pos=None, vec=None, target=False): - """ - Creates an Arc instance and transforms it to the desired position and - orientation. +def get_arc( + dogleg, radius, toolface, pos=None, vec=None, target=False +) -> tuple: + """Creates an Arc instance and transforms it to the desired position + and orientation. Parameters ---------- diff --git a/welleng/version.py b/welleng/version.py index 4ca39e7..b4e3540 100644 --- a/welleng/version.py +++ b/welleng/version.py @@ -1 +1 @@ -__version__ = '0.8.2' +__version__ = '0.8.3'