diff --git a/tests/test_iacs_rec34.py b/tests/test_iacs_rec34.py new file mode 100644 index 0000000..d6efbce --- /dev/null +++ b/tests/test_iacs_rec34.py @@ -0,0 +1,501 @@ +""" +Tests de Verificacion y Validacion segun IACS Recommendation No.34. + +IACS Rec.34, Rev.1 (2014) — "Verification and Validation of Marine Computer +Programs": + + par.3 Definitions: + - Verification: confirmar que el programa implementa correctamente + la formulacion matematica (comparacion con solucion analitica exacta). + - Validation: confirmar que la formulacion matematica describe + adecuadamente el fenomeno fisico real. + + par.4 Verification Methods: + par.4.3 — Comparacion con solucion analitica o semi-analitica. + par.4.4 — Prueba de convergencia de malla (mesh convergence). + par.4.5 — Prueba de simetria. + + par.5 Validation Methods: + par.5.2 — Comparacion con resultados experimentales publicados. + (Pendiente — requiere datos experimentales en sprint futuros) + + par.6 Documentacion requerida: + par.6.1 — Titulo, proposito, modulo. + par.6.2 — Algoritmos y referencias bibliograficas. + par.6.3 — Casos de prueba y tolerancias. + par.6.4 — Resultados esperados vs. obtenidos. + +Este modulo cubre la VERIFICACION (par.4) para el Modulo 1: + - Geometria NURBS (BSplineCurve, LoftedSurface) + - Hidrostáticos del casco (Hull) + - Serialización (.arsd round-trip) + +Tolerancias de aceptacion: + - Integrales directas (V, Awp): error relativo < 0.5 % con 41 estaciones + - Momentos de primer orden (LCB, KB): error relativo < 1 % + - Momentos de segundo orden (IT, IL): error relativo < 2 % + - Coeficientes adimensionales (Cb, Cw, Cm): error absoluto < 0.005 + - Serialización round-trip: error absoluto < 1e-9 (doble precision) + +Casco de verificacion: Wigley hull analitico (Wigley 1934) + y(x,z) = (B/2)[1-(2xi/L)^2][1-(zeta/T)^2] + donde xi = x-L/2, zeta = z-T + +Soluciones analiticas usadas: + V = 4BLT/9 + Cb = 4/9 + Awp(T) = 2BL/3 + Cw(T) = 2/3 + KB = 5T/8 + IT(T) = (2/3)(B/2)^3 (L/2)(32/35) + LCB = L/2 (simetria) + +Referencia: + Wigley, W.C.S. (1934). A comparison of experiment and calculated wave + profiles and wave resistances for a form having parabolic waterlines. + Proc. Roy. Soc. London A, 144, 144-159. + +Autor: Alvaro Romero | Modulo 1 -- AR-ShipDesign +IACS Rec.34 par.4.3, par.6 +""" +from __future__ import annotations + +import json +import math +import tempfile +from pathlib import Path + +import numpy as np +import pytest + +from arshipdesign.core.hull import Hull +from arshipdesign.core.offsets import OffsetsTable +from arshipdesign.core.project import Project +from arshipdesign.geometry.nurbs_curve import BSplineCurve +from arshipdesign.geometry.nurbs_surface import LoftedSurface + + +# --------------------------------------------------------------------------- +# Parametros del casco Wigley de verificacion +# (alta resolucion para minimizar error de cuadratura) +# --------------------------------------------------------------------------- + +L = 10.0 # Lpp [m] +B = 1.5 # manga [m] +T = 0.75 # calado [m] +NS = 41 # estaciones +NW = 21 # lineas de agua + +# Soluciones analiticas exactas +V_EXACT = 4.0 * B * L * T / 9.0 +CB_EXACT = 4.0 / 9.0 +AWP_EXACT = 2.0 * B * L / 3.0 +CW_EXACT = 2.0 / 3.0 +KB_EXACT = 5.0 * T / 8.0 +LCB_EXACT = L / 2.0 # simetria +IT_EXACT = (2.0/3.0)*(B/2)**3*(L/2)*(32.0/35.0) + +# Tolerancias IACS Rec.34 par.6.3 +TOL_V = 0.005 # 0.5 % error relativo para integrales de volumen +TOL_AWP = 0.005 # 0.5 % para area del plano de flotacion +TOL_COEF = 0.005 # 0.005 error absoluto para coeficientes adimensionales +TOL_KB = 0.010 # 1 % para centroide vertical (momento primer orden) +TOL_LCB = 0.005 # 0.5 % para LCB +TOL_IT = 0.020 # 2 % para segundo momento (IT) + + +@pytest.fixture(scope="module") +def wigley() -> Hull: + """Casco Wigley de referencia — resolucion alta para verificacion.""" + return Hull.from_wigley(lpp=L, beam=B, draft=T, + n_stations=NS, n_waterlines=NW) + + +# =========================================================================== +# A. Verificacion IACS Rec.34 par.4.3 — solucion analitica +# =========================================================================== + +class TestIACSVerificationAnalytic: + """Comparacion con la solucion analitica exacta del casco Wigley. + + IACS Rec.34 par.4.3: "Los resultados calculados se comparan con la + solucion analitica o semi-analitica del mismo problema." + """ + + def test_V001_volume_of_displacement(self, wigley: Hull) -> None: + """V001 — Volumen de desplazamiento V = 4BLT/9. + + Metodo: regla de Simpson sobre 41 secciones. + Tolerancia: < 0.5 % (IACS Rec.34 par.6.3). + """ + V = wigley.volume_of_displacement() + err_rel = abs(V - V_EXACT) / V_EXACT + assert err_rel < TOL_V, ( + f"V001 FAIL: V={V:.6f} m3, analitico={V_EXACT:.6f} m3, " + f"error_rel={err_rel*100:.3f}% > {TOL_V*100:.1f}%" + ) + + def test_V002_block_coefficient(self, wigley: Hull) -> None: + """V002 — Coeficiente de bloque Cb = V/(L*B*T) = 4/9. + + Tolerancia: |Cb_num - Cb_exact| < 0.005. + """ + cb = wigley.block_coefficient() + err = abs(cb - CB_EXACT) + assert err < TOL_COEF, ( + f"V002 FAIL: Cb={cb:.6f}, analitico={CB_EXACT:.6f}, " + f"error={err:.6f} > {TOL_COEF}" + ) + + def test_V003_waterplane_area(self, wigley: Hull) -> None: + """V003 — Area del plano de flotacion Awp = 2BL/3. + + Tolerancia: < 0.5 % error relativo. + """ + awp = wigley.waterplane_area() + err_rel = abs(awp - AWP_EXACT) / AWP_EXACT + assert err_rel < TOL_AWP, ( + f"V003 FAIL: Awp={awp:.6f} m2, analitico={AWP_EXACT:.6f} m2, " + f"error_rel={err_rel*100:.3f}%" + ) + + def test_V004_waterplane_coefficient(self, wigley: Hull) -> None: + """V004 — Coeficiente de plano de flotacion Cw = Awp/(L*B) = 2/3. + + Tolerancia: error absoluto < 0.005. + """ + cw = wigley.waterplane_coefficient() + err = abs(cw - CW_EXACT) + assert err < TOL_COEF, ( + f"V004 FAIL: Cw={cw:.6f}, analitico={CW_EXACT:.6f}, " + f"error={err:.6f}" + ) + + def test_V005_lcb_symmetry(self, wigley: Hull) -> None: + """V005 — LCB en el punto medio (L/2) por simetria longitudinal. + + El casco Wigley es simetrico respecto a x = L/2 => LCB = L/2. + Tolerancia: < 0.5 % de L. + """ + lcb = wigley.lcb() + err_rel = abs(lcb - LCB_EXACT) / L + assert err_rel < TOL_LCB, ( + f"V005 FAIL: LCB={lcb:.6f} m, analitico={LCB_EXACT:.6f} m, " + f"error_rel={err_rel*100:.3f}%" + ) + + def test_V006_vcb_kb(self, wigley: Hull) -> None: + """V006 — Centro vertical de carena KB = 5T/8. + + Tolerancia: < 1 % de T. + """ + kb = wigley.vcb() + err_rel = abs(kb - KB_EXACT) / T + assert err_rel < TOL_KB, ( + f"V006 FAIL: KB={kb:.6f} m, analitico={KB_EXACT:.6f} m, " + f"error_rel={err_rel*100:.3f}%" + ) + + def test_V007_it_waterplane(self, wigley: Hull) -> None: + """V007 — Segundo momento de area IT = (2/3)(B/2)^3(L/2)(32/35). + + Tolerancia: < 2 % error relativo. + """ + it = wigley.it_waterplane() + err_rel = abs(it - IT_EXACT) / IT_EXACT + assert err_rel < TOL_IT, ( + f"V007 FAIL: IT={it:.6f} m4, analitico={IT_EXACT:.6f} m4, " + f"error_rel={err_rel*100:.3f}%" + ) + + def test_V008_tpc_consistent_with_awp(self, wigley: Hull) -> None: + """V008 — TPC = Awp * rho / 100000 (verificacion de consistencia). + + Tolerancia: error relativo < 0.5 %. + """ + tpc = wigley.tpc(rho=1025.0) + tpc_ref = AWP_EXACT * 1025.0 / 100_000.0 + err_rel = abs(tpc - tpc_ref) / tpc_ref + assert err_rel < TOL_AWP, ( + f"V008 FAIL: TPC={tpc:.6f}, referencia={tpc_ref:.6f}, " + f"error_rel={err_rel*100:.3f}%" + ) + + def test_V009_km_transverse_chain(self, wigley: Hull) -> None: + """V009 — KMT = KB + BM_T = KB + IT/V (verificacion de cadena). + + Tolerancia: error absoluto < 1e-6 m (precision numerica interna). + """ + kmt = wigley.km_transverse() + kb = wigley.vcb() + it = wigley.it_waterplane() + v = wigley.volume_of_displacement() + kmt_ref = kb + it / v + assert abs(kmt - kmt_ref) < 1e-6, ( + f"V009 FAIL: KMT={kmt:.8f}, referencia={kmt_ref:.8f}" + ) + + +# =========================================================================== +# B. Verificacion IACS Rec.34 par.4.4 — convergencia de malla +# =========================================================================== + +class TestIACSMeshConvergence: + """Prueba de convergencia al refinar la discretizacion. + + IACS Rec.34 par.4.4: "Demostrar que los resultados convergen a un + valor estable al incrementar el numero de elementos." + """ + + @pytest.mark.parametrize("n_sta,n_wl", [ + (11, 6), + (21, 11), + (41, 21), + (81, 41), + ]) + def test_V010_volume_convergence(self, n_sta: int, n_wl: int) -> None: + """V010 — El volumen converge a V_EXACT al refinar la malla. + + La tolerancia se relaja para mallas gruesas: + n=11: < 5 % n=21: < 2 % n=41: < 0.5 % n=81: < 0.1 % + """ + hull = Hull.from_wigley(lpp=L, beam=B, draft=T, + n_stations=n_sta, n_waterlines=n_wl) + V = hull.volume_of_displacement() + err_rel = abs(V - V_EXACT) / V_EXACT + + tol_map = {11: 0.05, 21: 0.02, 41: 0.005, 81: 0.001} + tol = tol_map.get(n_sta, 0.05) + assert err_rel < tol, ( + f"V010 FAIL n={n_sta}: V={V:.5f} m3, error_rel={err_rel*100:.3f}% > {tol*100:.1f}%" + ) + + @pytest.mark.parametrize("n_sta,n_wl", [(11,6), (21,11), (41,21)]) + def test_V011_awp_convergence(self, n_sta: int, n_wl: int) -> None: + """V011 — Awp converge a 2BL/3 al refinar la malla.""" + hull = Hull.from_wigley(lpp=L, beam=B, draft=T, + n_stations=n_sta, n_waterlines=n_wl) + awp = hull.waterplane_area() + err_rel = abs(awp - AWP_EXACT) / AWP_EXACT + + tol_map = {11: 0.05, 21: 0.02, 41: 0.005} + tol = tol_map.get(n_sta, 0.05) + assert err_rel < tol, ( + f"V011 FAIL n={n_sta}: Awp={awp:.5f} m2, error_rel={err_rel*100:.3f}%" + ) + + def test_V012_volume_monotone_convergence(self) -> None: + """V012 — El error de volumen decrece al aumentar n (convergencia monotona). + + Verifica que el refinamiento siempre mejora el resultado. + """ + n_list = [11, 21, 41] + errors = [] + for n in n_list: + h = Hull.from_wigley(lpp=L, beam=B, draft=T, + n_stations=n, n_waterlines=n//2 + 1) + err = abs(h.volume_of_displacement() - V_EXACT) / V_EXACT + errors.append(err) + + for i in range(len(errors) - 1): + assert errors[i+1] <= errors[i] * 1.5, ( + f"V012 FAIL: error no decrece entre n={n_list[i]} y n={n_list[i+1]}: " + f"{errors[i]:.4e} -> {errors[i+1]:.4e}" + ) + + +# =========================================================================== +# C. Verificacion IACS Rec.34 par.4.5 — simetria +# =========================================================================== + +class TestIACSSymmetry: + """Prueba de simetria. + + IACS Rec.34 par.4.5: "Un problema simetrico debe producir una solucion + simetrica." + """ + + def test_V013_lcb_symmetry(self, wigley: Hull) -> None: + """V013 — LCB = L/2 para un casco simetrico longitudinalmente.""" + lcb = wigley.lcb() + assert abs(lcb - L / 2.0) / L < 1e-4, ( + f"V013 FAIL: LCB={lcb:.6f} m, esperado={L/2.0:.6f} m" + ) + + def test_V014_section_areas_symmetric(self, wigley: Hull) -> None: + """V014 — Las areas de cuadernas son simetricas respecto al midship. + + A(x) = A(L - x) para el casco Wigley. + """ + sections = wigley.offsets.to_sections() + n = len(sections) + for i in range(n // 2): + j = n - 1 - i + ai = sections[i].area(draft=T) + aj = sections[j].area(draft=T) + # Tolerancia 0.1 % (error de discretizacion) + if ai + aj > 1e-6: + err_rel = abs(ai - aj) / ((ai + aj) / 2.0) + assert err_rel < 0.001, ( + f"V014 FAIL i={i}: A[{i}]={ai:.5f}, A[{j}]={aj:.5f}, " + f"error_rel={err_rel*100:.4f}%" + ) + + def test_V015_offsets_table_symmetric(self, wigley: Hull) -> None: + """V015 — La tabla de offsets es simetrica respecto a x = L/2.""" + ot = wigley.offsets + n = ot.n_stations + for i in range(n // 2): + j = n - 1 - i + np.testing.assert_allclose( + ot.data[i, :], ot.data[j, :], + atol=1e-12, + err_msg=f"V015 FAIL: estaciones {i} y {j} no son simetricas" + ) + + +# =========================================================================== +# D. Verificacion de la geometria NURBS +# =========================================================================== + +class TestIACSNURBSVerification: + """Verificacion de los componentes de geometria NURBS. + + IACS Rec.34 par.4.3: comparacion con solucion analitica de la + longitud y posicion de puntos sobre curvas conocidas. + """ + + def test_V016_bspline_circle_approximation(self) -> None: + """V016 — BSplineCurve interpola exactamente puntos de una circunferencia.""" + # 9 puntos en un semicirculo de radio 1 + angles = np.linspace(0, math.pi, 9) + pts = np.column_stack([np.cos(angles), np.sin(angles)]) + curve = BSplineCurve(pts, degree=3) + # Verificar que los puntos de interpolacion estan en el semicirculo + for t_val in np.linspace(0, 1, 9): + p = curve.evaluate(t_val) + r = math.hypot(p[0], p[1]) + # Tolerancia 1 % (aprox. con grado 3) + assert abs(r - 1.0) < 0.05, ( + f"V016 FAIL: r={r:.4f} en t={t_val:.3f}" + ) + + def test_V017_bspline_line_exact(self) -> None: + """V017 — BSplineCurve es exacta para puntos colineales (linea recta).""" + pts = np.column_stack([np.linspace(0, 10, 7), np.zeros(7)]) + curve = BSplineCurve(pts, degree=3) + for t_val in np.linspace(0, 1, 20): + p = curve.evaluate(t_val) + # y debe ser 0, x debe estar en [0, 10] + assert abs(p[1]) < 1e-10, f"V017 FAIL: y={p[1]:.2e} en t={t_val:.3f}" + assert 0.0 - 1e-9 <= p[0] <= 10.0 + 1e-9 + + def test_V018_lofted_surface_wigley_midship(self) -> None: + """V018 — La superficie NURBS del Wigley es correcta en el midship. + + En x = L/2: la cuaderna debe tener semi-manga = B/2 * (1-(2*0/L)^2) + * (1-(zeta/T)^2) = (B/2) * f_zeta(z). + """ + hull = Hull.from_wigley(lpp=L, beam=B, draft=T, + n_stations=NS, n_waterlines=NW) + # Interpolacion bilineal en el punto medio + y_mid_top = hull.offsets.half_breadth(L / 2.0, T) + # En (x=L/2, z=T): f_xi = 1, f_zeta = 1 -> y = B/2 + assert abs(y_mid_top - B / 2.0) < 1e-9, ( + f"V018 FAIL: y(L/2,T)={y_mid_top:.6f}, esperado={B/2.0:.6f}" + ) + + def test_V019_lofted_surface_endpoints_zero(self) -> None: + """V019 — Semi-manga cero en AP y FP para el casco Wigley.""" + hull = Hull.from_wigley(lpp=L, beam=B, draft=T, + n_stations=NS, n_waterlines=NW) + # En x=0 (AP) y x=L (FP), la semi-manga debe ser 0 en todos los z + np.testing.assert_allclose(hull.offsets.data[0, :], 0.0, atol=1e-12) + np.testing.assert_allclose(hull.offsets.data[-1, :], 0.0, atol=1e-12) + + +# =========================================================================== +# E. Verificacion de la serializacion (trazabilidad de datos) +# =========================================================================== + +class TestIACSSerializationVerification: + """Verificacion de la serializacion segun IACS Rec.34 par.6.1. + + "Los datos de entrada (offsets) deben poder guardarse y restaurarse + sin perdida de precision." + """ + + def test_V020_serialization_preserves_volume(self, wigley: Hull) -> None: + """V020 — La serializacion preserva el volumen con precision doble.""" + d = wigley.to_dict() + h2 = Hull.from_dict(d) + assert abs(h2.volume_of_displacement() - + wigley.volume_of_displacement()) < 1e-9 + + def test_V021_project_roundtrip_preserves_it(self, wigley: Hull) -> None: + """V021 — IT se preserva exactamente tras guardar/cargar proyecto.""" + proj = Project.new("IACS V021") + proj.set_hull(wigley) + + with tempfile.TemporaryDirectory() as tmp: + p = Path(tmp) / "v021.arsd" + proj.save(p) + proj2 = Project.load(p) + h2 = proj2.hull + assert h2 is not None + assert abs(h2.it_waterplane() - wigley.it_waterplane()) < 1e-9 + + def test_V022_project_roundtrip_preserves_offsets_data( + self, wigley: Hull + ) -> None: + """V022 — La tabla de offsets es bit-a-bit identica tras round-trip.""" + proj = Project.new("IACS V022") + proj.set_hull(wigley) + with tempfile.TemporaryDirectory() as tmp: + p = Path(tmp) / "v022.arsd" + proj.save(p) + proj2 = Project.load(p) + h2 = proj2.hull + np.testing.assert_array_equal( + h2.offsets.data, wigley.offsets.data, + err_msg="V022 FAIL: tabla de offsets no identica tras round-trip" + ) + + def test_V023_json_human_readable(self, wigley: Hull) -> None: + """V023 — El JSON dentro del .arsd es legible (no binario, no base64). + + IACS Rec.34 par.6.1: "La documentacion debe permitir trazabilidad." + Los datos numericos deben ser legibles por un auditor externo. + """ + import zipfile + proj = Project.new("IACS V023") + proj.set_hull(wigley) + with tempfile.TemporaryDirectory() as tmp: + p = Path(tmp) / "v023.arsd" + proj.save(p) + with zipfile.ZipFile(p) as zf: + ship_txt = zf.read("ship.json").decode("utf-8") + # El texto debe contener floats legibles, no base64 + assert "hull_v1" in ship_txt + assert "x_stations" in ship_txt + # Debe ser JSON valido + ship = json.loads(ship_txt) + assert isinstance(ship["hull"]["offsets"]["x_stations"], list) + assert isinstance(ship["hull"]["offsets"]["x_stations"][0], float) + + +# =========================================================================== +# F. Resumen de cobertura IACS Rec.34 +# =========================================================================== + +class TestIACSCoverageSummary: + """Meta-test: verifica que todos los IDs de test V001-V023 existen.""" + + def test_all_verification_ids_covered(self) -> None: + """Comprueba que los IDs V001-V023 estan documentados en este modulo.""" + import inspect, sys + module = sys.modules[__name__] + src = inspect.getsource(module) + for i in range(1, 24): + vid = f"V{i:03d}" + assert vid in src, f"ID de verificacion '{vid}' no encontrado en el modulo"