From 3867ccfd706400bf6584e5ac4e669acf66f29c43 Mon Sep 17 00:00:00 2001 From: Matt Hall <5151457+mattEhall@users.noreply.github.com> Date: Wed, 20 Nov 2024 15:27:50 -0700 Subject: [PATCH 1/5] Preliminary support for lines out of the water: - Modifications in Line to apply reduced weight per unit length for Line objects with one end out of the water. These are a rough job and may best be replaced with an improved depth-aware Catenary function in the future. --- moorpy/Catenary.py | 6 +++++- moorpy/line.py | 27 +++++++++++++++++++++------ 2 files changed, 26 insertions(+), 7 deletions(-) diff --git a/moorpy/Catenary.py b/moorpy/Catenary.py index ed5d027..5db1053 100644 --- a/moorpy/Catenary.py +++ b/moorpy/Catenary.py @@ -1467,7 +1467,11 @@ def step_func_cat(X, args, Y, info, Ytarget, err, tols, iter, maxIter): #(fAH1, fAV1, fBH1, fBV1, info1) = catenary(0.16525209064463553, 1000.1008645356192, 976.86359774357, 861890783.955385, -5.663702119364548, 0.0, -0.0, 0.0, 17325264.3383085, 5.000000000000001e-05, 101, 1, 1000) #(fAH1, fAV1, fBH1, fBV1, info1) = catenary(0.0, 997.8909672113768, 1006.0699998926483, 361256000.00000006, 2.4952615384615333, 0, 0.0, 0, 0, 0.0001, 101, 100, 1, 1000) #(fAH1, fAV1, fBH1, fBV1, info1) = catenary(148.60942473375343, 1315.624259105325, 1279.7911844766932, 17497567581.802254, -81.8284437736437, CB=0.0, alpha=-0.0, HF0=795840.9590915733, VF0=6399974.444658389, Tol=0.0005, MaxIter=100, depth=1300.0, plots=1) - (fAH1, fAV1, fBH1, fBV1, info1) = catenary(2887.113193120885, -838.1577017360617, 2979.7592390255295, 41074520109.87981, 1104.878355564403, CB=0.0, alpha=81.41243787249468, HF0=0.0, VF0=3564574.0622291695, Tol=0.0005, MaxIter=100, depth=3000.0, plots=1) + #(fAH1, fAV1, fBH1, fBV1, info1) = catenary(2887.113193120885, -838.1577017360617, 2979.7592390255295, 41074520109.87981, 1104.878355564403, CB=0.0, alpha=81.41243787249468, HF0=0.0, VF0=3564574.0622291695, Tol=0.0005, MaxIter=100, depth=3000.0, plots=1) + + + (fAH1, fAV1, fBH1, fBV1, info1) = catenary(20.386176076554875, -83.33919733899663, 85.0, 469622033.24936056, -166.5252156894419, + CB=-897.8894086817654, alpha=0, HF0=1045436.0319202082, VF0=4280849.602990574, Tol=0.0001, MaxIter=100, depth=800, plots=1) #(fAH1, fAV1, fBH1, fBV1, info1) = catenary(274.9, 15.6, 328.5, 528887323., -121., # CB=-48., alpha=0, HF0=17939., VF0=20596., Tol=2e-05, MaxIter=100, plots=1) diff --git a/moorpy/line.py b/moorpy/line.py index ab285f0..ed352a7 100644 --- a/moorpy/line.py +++ b/moorpy/line.py @@ -743,6 +743,21 @@ def staticSolve(self, reset=False, tol=0.0001, profiles=0): else: self.cb = -depthA - self.rA[2] # when cb < 0, -cb is defined as height of end A off seabed (in catenary) + # Handle case of line above the water + if self.rA[2] > 0 or self.rB[2] > 0: # end B out of water << or do this in catenary? + ww = (self.type['m'] - np.pi/4*self.type['d_vol']**2 *self.sys.rho) * self.sys.g + wa = (self.type['m']) * self.sys.g # weight in air + + z_top = max(self.rA[2], self.rB[2]) + z_bot = min(self.rA[2], self.rB[2]) + + # overwrite the weight based on the wet-dry weight ratio based on end z vals + self.w = ww + (wa - ww) * ( -min(0, z_bot) )/( z_top - z_bot ) + else: + # reset to the normal weight if the line hasn't left the water + self.w = (self.type['m'] - np.pi/4*self.type['d_vol']**2 *self.sys.rho) * self.sys.g + + # >>> TODO: should replace all this with a more versatile catenary solve in future <<< # ----- Perform rotation/transformation to 2D plane of catenary ----- @@ -752,9 +767,9 @@ def staticSolve(self, reset=False, tol=0.0001, profiles=0): if np.sum(np.abs(self.fCurrent)) > 0: # total line exernal force per unit length vector (weight plus current drag) - w_vec = self.fCurrent/self.L + np.array([0, 0, -self.type["w"]]) - w = np.linalg.norm(w_vec) - w_hat = w_vec/w + w_vec = self.fCurrent/self.L + np.array([0, 0, -self.w]) + w_total = np.linalg.norm(w_vec) + w_hat = w_vec/w_total # >>> may need to adjust to handle case of buoyant vertical lines <<< @@ -773,7 +788,7 @@ def staticSolve(self, reset=False, tol=0.0001, profiles=0): # if no current force, things are simple else: R_curr = np.eye(3,3) - w = self.type["w"] + w_total = self.w # apply a rotation about Z' to align the line profile with the X'-Z' plane @@ -811,14 +826,14 @@ def staticSolve(self, reset=False, tol=0.0001, profiles=0): if 'EA' in self.type: try: (fAH, fAV, fBH, fBV, info) = catenary(LH, LV, self.L, self.EA, - w, CB=cb, alpha=alpha, HF0=self.HF, VF0=self.VF, Tol=tol, + w_total, CB=cb, alpha=alpha, HF0=self.HF, VF0=self.VF, Tol=tol, nNodes=self.nNodes, plots=profiles, depth=self.sys.depth) except CatenaryError as error: raise LineError(self.number, error.message) #If EA isnt found then we will use the ten-str relationship defined in the input file else: - (fAH, fAV, fBH, fBV, info) = nonlinear(LH, LV, self.L, self.type['Str'], self.type['Ten'],np.linalg.norm(w)) + (fAH, fAV, fBH, fBV, info) = nonlinear(LH, LV, self.L, self.type['Str'], self.type['Ten'],np.linalg.norm(w_total)) # save line profile coordinates in global frame (involves inverse rotation) From 3b1b85be421a5a7913c4012c331031475751ff3d Mon Sep 17 00:00:00 2001 From: Matt Hall <5151457+mattEhall@users.noreply.github.com> Date: Wed, 20 Nov 2024 15:37:33 -0700 Subject: [PATCH 2/5] More robust seabed slope support - Seabed slope for each Line is now calculated based on the depth beneath each endpoint (assuming a uniformly sloped seabed underneath the line) rather than looking at the panel slope under end A. - This should be more robust by avoiding inconsistencies when lines cross panels, though could reduce accuracy for finer grids with Lines that only slightly touch. --- moorpy/line.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/moorpy/line.py b/moorpy/line.py index ed352a7..30db83d 100644 --- a/moorpy/line.py +++ b/moorpy/line.py @@ -791,6 +791,11 @@ def staticSolve(self, reset=False, tol=0.0001, profiles=0): w_total = self.w + # horizontal and vertical dimensions of line profile (end A to B) + LH = np.linalg.norm(dr[:2]) + LV = dr[2] + + # apply a rotation about Z' to align the line profile with the X'-Z' plane theta_z = -np.arctan2(dr[1], dr[0]) R_z = rotationMatrix(0, 0, theta_z) @@ -802,8 +807,10 @@ def staticSolve(self, reset=False, tol=0.0001, profiles=0): if self.rA[2] <= -depthA or self.rB[2] <= -depthB: nvecA_prime = np.matmul(R, nvecA) - dz_dx = -nvecA_prime[0]*(1.0/nvecA_prime[2]) # seabed slope components - dz_dy = -nvecA_prime[1]*(1.0/nvecA_prime[2]) # seabed slope components + #dz_dx = -nvecA_prime[0]*(1.0/nvecA_prime[2]) # seabed slope components + #dz_dy = -nvecA_prime[1]*(1.0/nvecA_prime[2]) # seabed slope components + # Seabed slope along line direction (based on end A/B depth) + dz_dx = (-depthB + depthA)/LH # we only care about dz_dx since the line is in the X-Z plane in this rotated situation alpha = np.degrees(np.arctan(dz_dx)) cb = self.cb @@ -815,10 +822,6 @@ def staticSolve(self, reset=False, tol=0.0001, profiles=0): alpha = 0 cb = self.cb - # horizontal and vertical dimensions of line profile (end A to B) - LH = np.linalg.norm(dr[:2]) - LV = dr[2] - # ----- call catenary function or alternative and save results ----- From 969381db931d350f14ec6d597615125fb226b802 Mon Sep 17 00:00:00 2001 From: Matt Hall <5151457+mattEhall@users.noreply.github.com> Date: Fri, 22 Nov 2024 11:33:44 -0700 Subject: [PATCH 3/5] Support for lines laying fully along sloped seabed --- moorpy/Catenary.py | 175 ++++++++++++++++++++++++++++++--------------- 1 file changed, 116 insertions(+), 59 deletions(-) diff --git a/moorpy/Catenary.py b/moorpy/Catenary.py index 5db1053..eaee1ae 100644 --- a/moorpy/Catenary.py +++ b/moorpy/Catenary.py @@ -158,6 +158,7 @@ def catenary(XF, ZF, L, EA, W, CB=0, alpha=0, HF0=0, VF0=0, Tol=0.000001, # ProfileType=5: Similar to above but both ends are off seabed, so it's U shaped and fully slack # ProfileType=6: Completely vertical line without seabed contact (on the seabed is handled by 4 and 5) # ProfileType = 7: Portion of the line is resting on the seabed, and the seabed has a slope + # ProfileType = 8: line resting on the seabed, and the seabed has a slope EA_W = EA/W @@ -181,68 +182,120 @@ def dV_dZ_s(z0, H): # height off seabed to evaluate at (infinite if 0), horizo # because a large value here risks adding a bad cross coupling term in the system stiffness matrix # ProfileType 0 case - entirely along seabed - if ZF==0.0 and hA == 0.0 and W > 0: + if hB==0.0 and hA == 0.0 and W > 0: - # >>> this case should be extended to support sloped seabeds <<< - if not alpha == 0: raise CatenaryError("Line along seabed but seabed is sloped - not yet supported") - - ProfileType = 0 + if alpha == 0: # flat seabed case - if CB==0 or XF <= L: # case 1: no friction, or zero tension - HF = np.max([0, (XF/L - 1.0)*EA]) - HA = 1.0*HF - elif 0.5*L + EA/CB/W*(1-XF/L) <= 0: # case 2: seabed friction but tension at anchor (xB estimate < 0) - HF = (XF/L -1.0)*EA + 0.5*CB*W*L - HA = np.max([0.0, HF - CB*W*L]) - else: # case 3: seabed friction and zero anchor tension - if EA*CB*W*(XF-L) < 0: # something went wrong - breakpoint() - HF = np.sqrt(2*EA*CB*W*(XF-L)) - HA = 0.0 + ProfileType = 0 - VF = 0.0 - VA = 0.0 - - if HF > 0: # if taut - dHF_dXF = EA/L # approximation <<< what about friction? <<<<<<<< - #dVF_dZF = W + HF/L # vertical stiffness <<< approximation a - dVF_dZF = dV_dZ_s(Tol, HF) # vertical stiffness <<< approximation b - else: # if slack - dHF_dXF = 0.0 - dVF_dZF = W # vertical stiffness - - info["HF"] = HF # solution to be used to start next call (these are the solved variables, may be for anchor if line is reversed) - info["VF"] = 0.0 - info["stiffnessB"] = np.array([[ dHF_dXF, 0.0], [0.0, dVF_dZF]]) - info["stiffnessA"] = np.array([[ dHF_dXF, 0.0], [0.0, dVF_dZF]]) - info["stiffnessBA"] = np.array([[-dHF_dXF, 0.0], [0.0, 0.0]]) - info["LBot"] = L - info['ProfileType'] = 0 - info['Zextreme'] = 0 - - if plots > 0: + if CB==0 or XF <= L: # case 1: no friction, or zero tension + HF = np.max([0, (XF/L - 1.0)*EA]) + HA = 1.0*HF + elif 0.5*L + EA/CB/W*(1-XF/L) <= 0: # case 2: seabed friction but tension at anchor (xB estimate < 0) + HF = (XF/L -1.0)*EA + 0.5*CB*W*L + HA = np.max([0.0, HF - CB*W*L]) + else: # case 3: seabed friction and zero anchor tension + if EA*CB*W*(XF-L) < 0: # something went wrong + breakpoint() + HF = np.sqrt(2*EA*CB*W*(XF-L)) + HA = 0.0 + + VF = 0.0 + VA = 0.0 + + if HF > 0: # if taut + dHF_dXF = EA/L # approximation <<< what about friction? <<<<<<<< + #dVF_dZF = W + HF/L # vertical stiffness <<< approximation a + dVF_dZF = dV_dZ_s(Tol, HF) # vertical stiffness <<< approximation b + else: # if slack + dHF_dXF = 0.0 + dVF_dZF = W # vertical stiffness + + info["HF"] = HF # solution to be used to start next call (these are the solved variables, may be for anchor if line is reversed) + info["VF"] = 0.0 + info["stiffnessB"] = np.array([[ dHF_dXF, 0.0], [0.0, dVF_dZF]]) + info["stiffnessA"] = np.array([[ dHF_dXF, 0.0], [0.0, dVF_dZF]]) + info["stiffnessBA"] = np.array([[-dHF_dXF, 0.0], [0.0, 0.0]]) + info["LBot"] = L + info['ProfileType'] = 0 + info['Zextreme'] = 0 + + if plots > 0: + + if CB > 0 and XF > L: + xB = L - HF/W/CB # location of point at which line tension reaches zero + else: + xB = 0.0 + + # z values remain zero in this case + + if CB==0 or XF <= L: # case 1: no friction, or zero tension + Xs = XF/L*s # X values uniformly distributed + Te = Te + np.max([0, (XF/L - 1.0)*EA]) # uniform tension + elif xB <= 0: # case 2: seabed friction but tension at anchor + Xs = s*(1+CB*W/EA*(0.5*s-xB)) + Te = HF + CB*W*(s-L) + else: # case 3: seabed friction and zero anchor tension + for I in range(nNodes): + if s[I] <= xB: # if this node is in the zero tension range + Xs[I] = s[I]; # x is unstretched, z and Te remain zero + + else: # the tension is nonzero + Xs[I] = s[I] + CB*W/EA*(s[I] - xB)**2 + Te[I] = HF - CB*W*(L-s[I]) - if CB > 0 and XF > L: - xB = L - HF/W/CB # location of point at which line tension reaches zero - else: - xB = 0.0 + else: # sloped seabed case with full contact + + ProfileType = 8 - # z values remain zero in this case + cos_alpha = np.cos(np.radians(alpha)) + tan_alpha = np.tan(np.radians(alpha)) + sin_alpha = np.sin(np.radians(alpha)) - if CB==0 or XF <= L: # case 1: no friction, or zero tension - Xs = XF/L*s # X values uniformly distributed - Te = Te + np.max([0, (XF/L - 1.0)*EA]) # uniform tension - elif xB <= 0: # case 2: seabed friction but tension at anchor - Xs = s*(1+CB*W/EA*(0.5*s-xB)) - Te = HF + CB*W*(s-L) - else: # case 3: seabed friction and zero anchor tension - for I in range(nNodes): - if s[I] <= xB: # if this node is in the zero tension range - Xs[I] = s[I]; # x is unstretched, z and Te remain zero - - else: # the tension is nonzero - Xs[I] = s[I] + CB*W/EA*(s[I] - xB)**2 - Te[I] = HF - CB*W*(L-s[I]) + LAB = XF/cos_alpha # calling this the lenght along the seabed for this case + if LAB < 0: + breakpoint() + LAB = max(LAB, 0) # Ensure LAB is not less than zero + + + Ts = np.max([0, (LAB/L - 1.0)*EA]) # approximate tension based on streatch + Tg = L*W*sin_alpha # effect of weight and slope + # ignoring friction for now + + HF = cos_alpha*(Ts + 0.5*Tg) + VF = sin_alpha*(Ts + 0.5*Tg) + HA = cos_alpha*(Ts - 0.5*Tg) + VA = sin_alpha*(Ts - 0.5*Tg) + + if plots > 0: + + for I in range(nNodes): + Xs[I] = np.min([s[I]*cos_alpha, XF]) + Zs[I] = Xs[I]*tan_alpha + Te[I] = Ts - 0.5*Tg + Tg*float(I/(nNodes-1)) + + if HF > 0: # if taut + dHF_dXF = EA/L # CRUDE approximation + dVF_dZF = dV_dZ_s(Tol, HF) # same + else: # if slack + dHF_dXF = 0.0 + dVF_dZF = W # vertical stiffness + + info["HF"] = HF + info["VF"] = VF + info["stiffnessB"] = np.array([[ dHF_dXF, 0.0], [0.0, dVF_dZF]]) + info["stiffnessA"] = np.array([[ dHF_dXF, 0.0], [0.0, dVF_dZF]]) + info["stiffnessBA"] = np.array([[-dHF_dXF, 0.0], [0.0, 0.0]]) + info["LBot"] = L + info['ProfileType'] = 8 + info['Zextreme'] = 0 + + if plots > 0: + + if CB > 0 and XF > L: + xB = L - HF/W/CB # location of point at which line tension reaches zero + else: + xB = 0.0 # ProfileType 4 case - fully slack @@ -1187,7 +1240,7 @@ def eval_func_cat(X, args): info['error'] = True info['message'] = "ProfileType 7: VF_HF + SQRT1VF_HF2 <= 0" - elif HF*1000000 < VF: + elif HF*1000000000 < VF: info['error'] = True info['message'] = "ProfileType 7: HF << VF, line is slack, not supported yet" breakpoint() @@ -1470,8 +1523,12 @@ def step_func_cat(X, args, Y, info, Ytarget, err, tols, iter, maxIter): #(fAH1, fAV1, fBH1, fBV1, info1) = catenary(2887.113193120885, -838.1577017360617, 2979.7592390255295, 41074520109.87981, 1104.878355564403, CB=0.0, alpha=81.41243787249468, HF0=0.0, VF0=3564574.0622291695, Tol=0.0005, MaxIter=100, depth=3000.0, plots=1) - (fAH1, fAV1, fBH1, fBV1, info1) = catenary(20.386176076554875, -83.33919733899663, 85.0, 469622033.24936056, -166.5252156894419, - CB=-897.8894086817654, alpha=0, HF0=1045436.0319202082, VF0=4280849.602990574, Tol=0.0001, MaxIter=100, depth=800, plots=1) + #(fAH1, fAV1, fBH1, fBV1, info1) = catenary(20.386176076554875, -83.33919733899663, 85.0, 469622033.24936056, -166.5252156894419, + # CB=-897.8894086817654, alpha=0, HF0=1045436.0319202082, VF0=4280849.602990574, Tol=0.0001, MaxIter=100, depth=800, plots=1) + + + (fAH1, fAV1, fBH1, fBV1, info1) = catenary(242.43063215088046, 16.438515754164882, 242.94180148966598, 1977564385.9456, 3941.831324315529, CB=0.0, alpha=3.879122219600495, HF0=400446.3427158061, VF0=26452.686481212237, Tol=0.0001, MaxIter=100, depth=820.2155229016184, plots=1) + #(fAH1, fAV1, fBH1, fBV1, info1) = catenary(113.98763202381406, 135.64003587014417, 222.9088982113214, 1799620189.0374997, 3587.1279256290677, CB=0.0, alpha=-1.6131140432335964, HF0=737382.1987996304, VF0=1014826.8890259801, Tol=0.0001, MaxIter=100, depth=914.6515127916663, plots=1) #(fAH1, fAV1, fBH1, fBV1, info1) = catenary(274.9, 15.6, 328.5, 528887323., -121., # CB=-48., alpha=0, HF0=17939., VF0=20596., Tol=2e-05, MaxIter=100, plots=1) From 8973cfc025ce8e214ce93c55147fc5a7c36c7402 Mon Sep 17 00:00:00 2001 From: Matt Hall <5151457+mattEhall@users.noreply.github.com> Date: Mon, 23 Dec 2024 10:25:24 -0700 Subject: [PATCH 4/5] Line and Subsystem shadows now plot on seabed - Adjusted drawLine methdos to now draw the shadow on the seabed bathymetry grid if it exists. --- moorpy/line.py | 11 ++++++++++- moorpy/subsystem.py | 15 ++++++++++++--- 2 files changed, 22 insertions(+), 4 deletions(-) diff --git a/moorpy/line.py b/moorpy/line.py index 30db83d..df278fb 100644 --- a/moorpy/line.py +++ b/moorpy/line.py @@ -564,7 +564,16 @@ def drawLine(self, Time, ax, color="k", endpoints=False, shadow=True, colortensi linebit.append(ax.plot(Xs, Ys, Zs, color=color, lw=lw, zorder=100)) if shadow: - ax.plot(Xs, Ys, np.zeros_like(Xs)-self.sys.depth, color=[0.5, 0.5, 0.5, 0.2], lw=lw, zorder = 1.5) # draw shadow + if self.sys.seabedMod == 0: + Zs = np.zeros_like(Xs)-self.sys.depth + elif self.seabedMod == 1: + Zs = self.sys.depth - self.sys.xSlope*Xs - self.sys.ySlope*Ys + elif self.seabedMod == 2: + Zs = np.zeros(len(Xs)) + for i in range(len(Xs)): + Zs[i] = self.sys.getDepthAtLocation(Xs[i], Ys[i]) + + ax.plot(Xs, Ys, Zs, color=[0.5, 0.5, 0.5, 0.2], lw=lw, zorder = 1.5) # draw shadow if endpoints == True: linebit.append(ax.scatter([Xs[0], Xs[-1]], [Ys[0], Ys[-1]], [Zs[0], Zs[-1]], color = color)) diff --git a/moorpy/subsystem.py b/moorpy/subsystem.py index 782c748..e3a5e3e 100644 --- a/moorpy/subsystem.py +++ b/moorpy/subsystem.py @@ -284,7 +284,7 @@ def setEndPosition(self, r, endB, sink=False): raise LineError("setEndPosition: endB value has to be either 1 or 0") - def staticSolve(self, reset=False, tol=0, profiles=0, maxIter = 500): + def staticSolve(self, tol=0, profiles=0, maxIter = 500): '''Solve internal equilibrium of the Subsystem and saves the forces and stiffnesses at the ends in the global reference frame. All the This method mimics the behavior of the Line.staticSolve method and @@ -292,7 +292,7 @@ def staticSolve(self, reset=False, tol=0, profiles=0, maxIter = 500): equilibrium happens in the local 2D plane. Values in this local frame are also saved. ''' - +# >>> need to deal with bathymetry transformation for subsystems! >>> if tol==0: tol=self.eqtol @@ -462,7 +462,16 @@ def drawLine(self, Time, ax, color="k", endpoints=False, shadow=True, colortensi ax.plot(Xs, Ys, Zs, color=line.color, lw=line.lw, zorder=100) if shadow: - ax.plot(Xs, Ys, np.zeros_like(Xs)-self.depth, color=[0.5, 0.5, 0.5, 0.2], lw=line.lw, zorder = 1.5) # draw shadow + if self.seabedMod == 0: + Zs = np.zeros_like(Xs)-self.depth + elif self.seabedMod == 1: + Zs = self.depth - self.xSlope*Xs - self.ySlope*Ys + elif self.seabedMod == 2: + Zs = np.zeros(len(Xs)) + for i in range(len(Xs)): + Zs[i] = self.getDepthAtLocation(Xs[i], Ys[i]) + + ax.plot(Xs, Ys, Zs, color=[0.5, 0.5, 0.5, 0.2], lw=line.lw, zorder = 1.5) # draw shadow if endpoints == True: #linebit.append(ax.scatter([Xs[0], Xs[-1]], [Ys[0], Ys[-1]], [Zs[0], Zs[-1]], color = color)) From f55b264f581a2e37ae23667d2b50d8b15c0dc038 Mon Sep 17 00:00:00 2001 From: Stein Date: Fri, 22 Nov 2024 11:59:44 -0700 Subject: [PATCH 5/5] Subsystem bathymetry functionality - Included a function in subsystem.py to pull from the global bathymetry and produce a 2D bathymetry directly underneath the subsystem - - The bathGrid variable now becomes a 1D array of depths along the subsystem - - bathGrid_Ys become 0, and the xs are set as the nodes along each line in the subsystem - - This way, the alpha (slope) is calculated for each line in the subsystem's line list properly - This function is called every time subsystem.staticSolve is called, with updated bathymetry based on any new rA/rB positions - Additional work to do - - Make sure the global bathymetry is tracked properly (currently needs self.sys != None) - - Ensure the seabedMod parameter is set up properly --- moorpy/subsystem.py | 66 +++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 64 insertions(+), 2 deletions(-) diff --git a/moorpy/subsystem.py b/moorpy/subsystem.py index e3a5e3e..1064fbf 100644 --- a/moorpy/subsystem.py +++ b/moorpy/subsystem.py @@ -107,9 +107,21 @@ def __init__(self, mooringSys=None, num=0, depth=0, rho=1025, g=9.81, self.xSlope = getFromDict(kwargs, 'xSlope', default=0) self.ySlope = getFromDict(kwargs, 'ySlope', default=0) + #if 'bathymetry' in kwargs: + #self.seabedMod = 2 + #self.bathGrid_Xs, self.bathGrid_Ys, self.bathGrid = self.readBathymetryFile(kwargs['bathymetry']) + ''' if 'bathymetry' in kwargs: self.seabedMod = 2 - self.bathGrid_Xs, self.bathGrid_Ys, self.bathGrid = self.readBathymetryFile(kwargs['bathymetry']) + if isinstance(kwargs['bathymetry'],str): + self.bathGrid_Xs, self.bathGrid_Ys, self.bathGrid = readBathymetryFile(kwargs['bathymetry']) + elif isinstance(kwargs['bathymetry'],dict): + bath_dictionary = kwargs['bathymetry'] + # grid sent in, just assign to properties + self.bathGrid_Xs = bath_dictionary['x'] + self.bathGrid_Ys = bath_dictionary['y'] + self.bathGrid = bath_dictionary['depth'] + ''' # Note, System and Subsystem bathymetry can be set after creation # by setBathymetry, which uses link to existing data, for efficiency. @@ -124,6 +136,53 @@ def __init__(self, mooringSys=None, num=0, depth=0, rho=1025, g=9.81, self.MDoptions = {} # dictionary that can hold any MoorDyn options read in from an input file, so they can be saved in a new MD file if need be self.qs = 1 # flag that it's a MoorPy analysis, so System methods don't complain. One day should replace this <<< + + def setSSBathymetry(self): + '''Provides global bathymetry data, but creates a 2D line corresponding to what's underneath the subsystem + + - not naming it 'setBathymetry' since that is already a function name in the System class, but maybe I could name it the same? + + - should theoretically be called every iteration / every call to a new position, to get the most update bathymetry plane + + - need to make sure that any bathymetric operation in subsystem uses this new bathymetry information, and not the System.getDepthFromBathymetry + - - the fact that the bathGrid_Xs and Ys and grid arrays are new, should solve this issue + + - could set up something for seabedMod = 0, 1, 2 + + Parameters + ---------- + x : list or array + X coordinates of a global, rectangular bathymetry grid [m]. + y : list or array + Y coordinates of a global, rectangular bathymetry grid [m]. + depth : 2D array + Matrix of depths at x and y locations [m]. Positive depths. + Shape should be [len(x), len(y)] + ''' + + # get the depth at the anchor point of the subsystem + depthA, nvecA = self.sys.getDepthFromBathymetry(self.rA[0], self.rA[1]) + # get the depth under the fairlead point of the subsystem + depthB, nvecB = self.sys.getDepthFromBathymetry(self.rB[0], self.rB[1]) + + # collection information to create a new list of x coordinates to use for bathymetry for the subsystem (subsystems don't have y coordinates) + nNodes = [line.nNodes for line in self.lineList] # collecting the number of nodes in all the line sections + nx = np.sum(nNodes) - len(self.lineList) + 1 # adding all those nodes, subtracting the top nodes from each line (to avoid duplicates) and adding the fairlead node + xs = np.linspace(-self.span, 0, nx) # assuming that the subsystem rB will always be at 0,0 + + # calculate the bathymetric depths directly underneath the subsystem + depths = np.zeros(len(xs)) + for i in range(len(xs)): + x = self.rA[0] + (xs[i]+self.span)*self.cos_th # take the global x position of the anchor and add the cosine of the node x (while setting it relative to 0,0) + y = self.rA[1] + (xs[i]+self.span)*self.sin_th # take the global y position of the anchor and add the sine of the node x (while setting it relative to 0,0) + depths[i], _ = self.sys.getDepthFromBathymetry(x, y) # calculate the depth at this individual point for one row of depths + + # save all these variables to use later (overwrites any previous 'System' bathymetry arrays, and now getDepthFromBathymetry can use this data) + self.bathGrid_Xs = np.array(xs) + self.bathGrid_Ys = np.array([0]) + self.bathGrid = np.array([depths]) + + self.seabedMod = 2 def initialize(self, daf_dict={}): @@ -284,7 +343,7 @@ def setEndPosition(self, r, endB, sink=False): raise LineError("setEndPosition: endB value has to be either 1 or 0") - def staticSolve(self, tol=0, profiles=0, maxIter = 500): + def staticSolve(self, tol=0, profiles=0, maxIter=500): '''Solve internal equilibrium of the Subsystem and saves the forces and stiffnesses at the ends in the global reference frame. All the This method mimics the behavior of the Line.staticSolve method and @@ -308,6 +367,9 @@ def staticSolve(self, tol=0, profiles=0, maxIter = 500): else: self.pointList[ 0].setPosition([ -self.span , 0, self.rA[2]]) self.pointList[-1].setPosition([ -self.span+LH, 0, self.rB[2]]) + + # update the bathymetry based on any new positions + #self.setSSBathymetry() # get equilibrium self.solveEquilibrium(tol=tol, maxIter = maxIter)