From 1164066abf27f8c9b6e6253eb7f0fe8599d78a9b Mon Sep 17 00:00:00 2001 From: Piotr Bartman Date: Wed, 3 May 2023 22:25:00 +0200 Subject: [PATCH 01/34] BBO: init --- examples/Bagirov_Bartman_Ochal_2023.py | 109 +++++++++++++++++++++++++ 1 file changed, 109 insertions(+) create mode 100644 examples/Bagirov_Bartman_Ochal_2023.py diff --git a/examples/Bagirov_Bartman_Ochal_2023.py b/examples/Bagirov_Bartman_Ochal_2023.py new file mode 100644 index 00000000..104f7cfa --- /dev/null +++ b/examples/Bagirov_Bartman_Ochal_2023.py @@ -0,0 +1,109 @@ +""" +Created at 21.08.2019 +""" +from dataclasses import dataclass + +import numpy as np +from conmech.helpers.config import Config +from conmech.mesh.boundaries_description import BoundariesDescription +from conmech.plotting.drawer import Drawer +from conmech.scenarios.problems import ContactLaw, Static +from conmech.simulations.problem_solver import Static as StaticProblemSolver + + +class JureczkaOchal2019(ContactLaw): + @staticmethod + def potential_normal_direction(u_nu: float) -> float: + # if u_nu <= 0: + # return 0.0 + # if u_nu < 0.1: + # return 10 * u_nu * u_nu + # return 0.1 + if u_nu >= 0: + return 0.0 + if u_nu > -0.5e-3: + return 30e3 * u_nu ** 2 + if u_nu > -1e-3: + return 10e3 * (u_nu ** 2 + u_nu) + if u_nu > -2e-3: + return 5e3 * (u_nu ** 2 + u_nu) + return 0.0 + + @staticmethod + def potential_tangential_direction(u_tau: np.ndarray) -> float: + return np.log(np.sum(u_tau * u_tau) ** 0.5 + 1) + + @staticmethod + def subderivative_normal_direction(u_nu: float, v_nu: float) -> float: + return 0 + + @staticmethod + def regularized_subderivative_tangential_direction( + u_tau: np.ndarray, v_tau: np.ndarray, rho=1e-7 + ) -> float: + """ + Coulomb regularization + """ + return 0 + # regularization = 1 / np.sqrt(u_tau[0] * u_tau[0] + u_tau[1] * u_tau[1] + rho ** 2) + # result = regularization * (u_tau[0] * v_tau[0] + u_tau[1] * v_tau[1]) + # return result + + +mesh_density = 4 + + +@dataclass() +class StaticSetup(Static): + grid_height: ... = 0.1 + elements_number: ... = (mesh_density, 8 * mesh_density) + mu_coef: ... = 5.3e4 + la_coef: ... = 7.95e4 + contact_law: ... = JureczkaOchal2019 + + @staticmethod + def inner_forces(x, t=None): + return np.array([0.0, 0.0]) + + @staticmethod + def outer_forces(x, t=None): + return np.array([0, 5]) + + @staticmethod + def friction_bound(u_nu: float) -> float: + # if u_nu < 0: + # return 0 + # if u_nu < 0.1: + # return 8 * u_nu + return 0.0 + + boundaries: ... = BoundariesDescription( + contact=lambda x: x[1] == 0, dirichlet=lambda x: x[0] == 0 + ) + + +def main(show: bool = True, save: bool = False): + setup = StaticSetup(mesh_type="cross") + + for force in np.arange(20, 30 + 1, 2): + def outer_forces(x, t=None): + return np.array([0, force]) + + setup.outer_forces = outer_forces + + runner = StaticProblemSolver(setup, "schur") + + state = runner.solve( + verbose=True, + fixed_point_abs_tol=0.001, + initial_displacement=setup.initial_displacement, + method="Powell" + ) + config = Config() + drawer = Drawer(state=state, config=config) + drawer.colorful = True + drawer.draw(show=show, save=save) + + +if __name__ == "__main__": + main(show=True, save=False) From 424e3218c3062a1ed96aba0b3513111af09370c1 Mon Sep 17 00:00:00 2001 From: Piotr Bartman Date: Wed, 10 May 2023 18:43:01 +0200 Subject: [PATCH 02/34] BBO: WIP --- conmech/plotting/drawer.py | 2 +- examples/Bagirov_Bartman_Ochal_2023.py | 77 +++++++++++++------------- 2 files changed, 40 insertions(+), 39 deletions(-) diff --git a/conmech/plotting/drawer.py b/conmech/plotting/drawer.py index 3d7fe7e5..e013325c 100644 --- a/conmech/plotting/drawer.py +++ b/conmech/plotting/drawer.py @@ -76,7 +76,7 @@ def draw( # plt.axis("on") # axes.tick_params(left=True, bottom=True, labelleft=True, labelbottom=True) # - # axes.set_aspect("equal", adjustable="box") + axes.set_aspect("equal", adjustable="box") if title is not None: plt.title(title) diff --git a/examples/Bagirov_Bartman_Ochal_2023.py b/examples/Bagirov_Bartman_Ochal_2023.py index 104f7cfa..3abe63b4 100644 --- a/examples/Bagirov_Bartman_Ochal_2023.py +++ b/examples/Bagirov_Bartman_Ochal_2023.py @@ -11,23 +11,27 @@ from conmech.simulations.problem_solver import Static as StaticProblemSolver -class JureczkaOchal2019(ContactLaw): +class MMLV99(ContactLaw): @staticmethod def potential_normal_direction(u_nu: float) -> float: - # if u_nu <= 0: + # if u_nu >= 0: # return 0.0 - # if u_nu < 0.1: - # return 10 * u_nu * u_nu - # return 0.1 + # if u_nu > -0.5e-3: + # return (30 / 200) * u_nu ** 2 + # if u_nu > -1e-3: + # return (10 / 200) * (u_nu ** 2 + u_nu) + # if u_nu > -2e-3: + # return (5 / 200) * (u_nu ** 2 + u_nu + 2) + # return (40 / 200) if u_nu >= 0: return 0.0 - if u_nu > -0.5e-3: - return 30e3 * u_nu ** 2 + if u_nu > -0.5e-6: + return (30e6) * u_nu ** 2 if u_nu > -1e-3: - return 10e3 * (u_nu ** 2 + u_nu) + return (10e6) * (u_nu ** 2 + u_nu) - 4995 if u_nu > -2e-3: - return 5e3 * (u_nu ** 2 + u_nu) - return 0.0 + return (5e6) * (u_nu ** 2 + u_nu) + 10 + return 10030 @staticmethod def potential_tangential_direction(u_tau: np.ndarray) -> float: @@ -45,9 +49,6 @@ def regularized_subderivative_tangential_direction( Coulomb regularization """ return 0 - # regularization = 1 / np.sqrt(u_tau[0] * u_tau[0] + u_tau[1] * u_tau[1] + rho ** 2) - # result = regularization * (u_tau[0] * v_tau[0] + u_tau[1] * v_tau[1]) - # return result mesh_density = 4 @@ -55,11 +56,11 @@ def regularized_subderivative_tangential_direction( @dataclass() class StaticSetup(Static): - grid_height: ... = 0.1 + grid_height: ... = 0.01 elements_number: ... = (mesh_density, 8 * mesh_density) - mu_coef: ... = 5.3e4 - la_coef: ... = 7.95e4 - contact_law: ... = JureczkaOchal2019 + mu_coef: ... = (1.378e8) / (2 * (1 + 0.3)) + la_coef: ... = ((1.378e8) * 0.3) / ((1 + 0.3) * (1 - 2 * 0.3)) + contact_law: ... = MMLV99 @staticmethod def inner_forces(x, t=None): @@ -71,10 +72,6 @@ def outer_forces(x, t=None): @staticmethod def friction_bound(u_nu: float) -> float: - # if u_nu < 0: - # return 0 - # if u_nu < 0.1: - # return 8 * u_nu return 0.0 boundaries: ... = BoundariesDescription( @@ -82,28 +79,32 @@ def friction_bound(u_nu: float) -> float: ) -def main(show: bool = True, save: bool = False): +def main(save: bool = False): setup = StaticSetup(mesh_type="cross") - for force in np.arange(20, 30 + 1, 2): - def outer_forces(x, t=None): - return np.array([0, force]) + for method in ("Powell", "BFGS", "qsm")[2:]: + for force in np.arange(80000, 200000 + 1, 10000): + def outer_forces(x, t=None): + if x[1] >= 0.0099: + return np.array([0, force]) + return np.array([0, 0]) + - setup.outer_forces = outer_forces + setup.outer_forces = outer_forces - runner = StaticProblemSolver(setup, "schur") + runner = StaticProblemSolver(setup, "schur") - state = runner.solve( - verbose=True, - fixed_point_abs_tol=0.001, - initial_displacement=setup.initial_displacement, - method="Powell" - ) - config = Config() - drawer = Drawer(state=state, config=config) - drawer.colorful = True - drawer.draw(show=show, save=save) + state = runner.solve( + verbose=True, + fixed_point_abs_tol=0.001, + initial_displacement=setup.initial_displacement, + method=method + ) + config = Config() + drawer = Drawer(state=state, config=config) + drawer.colorful = True + drawer.draw(show=not save, save=save, title=f"{method}: {force}") if __name__ == "__main__": - main(show=True, save=False) + main(save=True) From ebaac665611ff440b37b02b9ab26da877a5683b3 Mon Sep 17 00:00:00 2001 From: Piotr Bartman Date: Thu, 11 May 2023 00:53:16 +0200 Subject: [PATCH 03/34] BBO: WIP x 8 --- examples/Bagirov_Bartman_Ochal_2023.py | 38 +++++++++++++++++--------- 1 file changed, 25 insertions(+), 13 deletions(-) diff --git a/examples/Bagirov_Bartman_Ochal_2023.py b/examples/Bagirov_Bartman_Ochal_2023.py index 3abe63b4..af5d600b 100644 --- a/examples/Bagirov_Bartman_Ochal_2023.py +++ b/examples/Bagirov_Bartman_Ochal_2023.py @@ -23,15 +23,16 @@ def potential_normal_direction(u_nu: float) -> float: # if u_nu > -2e-3: # return (5 / 200) * (u_nu ** 2 + u_nu + 2) # return (40 / 200) - if u_nu >= 0: + u_nu = -u_nu + if u_nu <= 0: return 0.0 - if u_nu > -0.5e-6: - return (30e6) * u_nu ** 2 - if u_nu > -1e-3: - return (10e6) * (u_nu ** 2 + u_nu) - 4995 - if u_nu > -2e-3: - return (5e6) * (u_nu ** 2 + u_nu) + 10 - return 10030 + if u_nu < 0.5 * mm: + return (30e6 * kN * surface) * u_nu ** 2 + if u_nu < 1 * mm: + return (10e6 * kN * surface) * (u_nu ** 2 + u_nu) - 4995 * surface * 1000 + if u_nu < 2 * mm: + return (5e6 * kN * surface) * (u_nu ** 2 + u_nu) + 10 * surface * 1000 + return 10030 * surface * 1000 @staticmethod def potential_tangential_direction(u_tau: np.ndarray) -> float: @@ -52,14 +53,19 @@ def regularized_subderivative_tangential_direction( mesh_density = 4 +kN = 1000 +mm = 0.001 +E = 1.378e8 * kN +kappa = 0.3 +surface = 5 * mm * 80 * mm * 8 @dataclass() class StaticSetup(Static): - grid_height: ... = 0.01 + grid_height: ... = 10 * mm elements_number: ... = (mesh_density, 8 * mesh_density) - mu_coef: ... = (1.378e8) / (2 * (1 + 0.3)) - la_coef: ... = ((1.378e8) * 0.3) / ((1 + 0.3) * (1 - 2 * 0.3)) + mu_coef: ... = (E) / (2 * (1 + kappa)) + la_coef: ... = ((E) * kappa) / ((1 + kappa) * (1 - 2 * kappa)) contact_law: ... = MMLV99 @staticmethod @@ -83,13 +89,12 @@ def main(save: bool = False): setup = StaticSetup(mesh_type="cross") for method in ("Powell", "BFGS", "qsm")[2:]: - for force in np.arange(80000, 200000 + 1, 10000): + for force in np.arange(20e3 * kN, 30e3 * kN + 1, 1e3 * kN): def outer_forces(x, t=None): if x[1] >= 0.0099: return np.array([0, force]) return np.array([0, 0]) - setup.outer_forces = outer_forces runner = StaticProblemSolver(setup, "schur") @@ -107,4 +112,11 @@ def outer_forces(x, t=None): if __name__ == "__main__": + from matplotlib import pyplot as plt + X = np.linspace(0, -3 * mm, 1000) + Y = np.empty(1000) + for i in range(1000): + Y[i] = MMLV99.potential_normal_direction(X[i]) + plt.plot(X, Y) + plt.show() main(save=True) From 59e227efcbdca7ef507c022a6811482ea2fa31ac Mon Sep 17 00:00:00 2001 From: Piotr Bartman Date: Thu, 11 May 2023 10:38:29 +0200 Subject: [PATCH 04/34] BBO: reproduced --- examples/Bagirov_Bartman_Ochal_2023.py | 55 +++++++++++++++----------- 1 file changed, 31 insertions(+), 24 deletions(-) diff --git a/examples/Bagirov_Bartman_Ochal_2023.py b/examples/Bagirov_Bartman_Ochal_2023.py index af5d600b..a4879346 100644 --- a/examples/Bagirov_Bartman_Ochal_2023.py +++ b/examples/Bagirov_Bartman_Ochal_2023.py @@ -9,30 +9,45 @@ from conmech.plotting.drawer import Drawer from conmech.scenarios.problems import ContactLaw, Static from conmech.simulations.problem_solver import Static as StaticProblemSolver +mesh_density = 4 +kN = 1000 +mm = 0.001 +E = 1.378e8 * kN +kappa = 0.3 +surface = 5 * mm * 80 * mm +k0 = (30e6 * kN * surface) +k10 = (10e6 * kN * surface) +k11 = (10e3 * kN * surface) +k20 = (5e6 * kN * surface) +k21 = (5e3 * kN * surface) + + +def normal_direction(u_nu: float) -> float: + u_nu = -u_nu + if u_nu <= 0: + return 0.0 + if u_nu < 0.5 * mm: + return k0 * u_nu * 2 + if u_nu < 1 * mm: + return k10 * (u_nu * 2) + k11 + if u_nu < 2 * mm: + return k20 * (u_nu * 2) + k21 + return 0 class MMLV99(ContactLaw): @staticmethod def potential_normal_direction(u_nu: float) -> float: - # if u_nu >= 0: - # return 0.0 - # if u_nu > -0.5e-3: - # return (30 / 200) * u_nu ** 2 - # if u_nu > -1e-3: - # return (10 / 200) * (u_nu ** 2 + u_nu) - # if u_nu > -2e-3: - # return (5 / 200) * (u_nu ** 2 + u_nu + 2) - # return (40 / 200) u_nu = -u_nu if u_nu <= 0: return 0.0 if u_nu < 0.5 * mm: - return (30e6 * kN * surface) * u_nu ** 2 + return k0 * u_nu ** 2 if u_nu < 1 * mm: - return (10e6 * kN * surface) * (u_nu ** 2 + u_nu) - 4995 * surface * 1000 + return k10 * u_nu ** 2 + k11 * u_nu if u_nu < 2 * mm: - return (5e6 * kN * surface) * (u_nu ** 2 + u_nu) + 10 * surface * 1000 - return 10030 * surface * 1000 + return k20 * u_nu ** 2 + k21 * u_nu + 4.0 + return 16 @staticmethod def potential_tangential_direction(u_tau: np.ndarray) -> float: @@ -52,20 +67,12 @@ def regularized_subderivative_tangential_direction( return 0 -mesh_density = 4 -kN = 1000 -mm = 0.001 -E = 1.378e8 * kN -kappa = 0.3 -surface = 5 * mm * 80 * mm * 8 - - @dataclass() class StaticSetup(Static): grid_height: ... = 10 * mm elements_number: ... = (mesh_density, 8 * mesh_density) - mu_coef: ... = (E) / (2 * (1 + kappa)) - la_coef: ... = ((E) * kappa) / ((1 + kappa) * (1 - 2 * kappa)) + mu_coef: ... = (E * surface) / (2 * (1 + kappa)) + la_coef: ... = ((E * surface) * kappa) / ((1 + kappa) * (1 - 2 * kappa)) contact_law: ... = MMLV99 @staticmethod @@ -89,7 +96,7 @@ def main(save: bool = False): setup = StaticSetup(mesh_type="cross") for method in ("Powell", "BFGS", "qsm")[2:]: - for force in np.arange(20e3 * kN, 30e3 * kN + 1, 1e3 * kN): + for force in np.arange(25e3 * kN, 26e3 * kN + 1, 1e3 * kN) * surface: def outer_forces(x, t=None): if x[1] >= 0.0099: return np.array([0, force]) From 3cfa9d88f53073f30cbaaa619e99d76377e1204d Mon Sep 17 00:00:00 2001 From: Piotr Bartman Date: Thu, 11 May 2023 20:04:35 +0200 Subject: [PATCH 05/34] BBO: reproduced --- examples/Bagirov_Bartman_Ochal_2023.py | 4 ++-- examples/Sofonea_Ochal_Bartman_2023.py | 6 ++++++ 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/examples/Bagirov_Bartman_Ochal_2023.py b/examples/Bagirov_Bartman_Ochal_2023.py index a4879346..d3e2ae2e 100644 --- a/examples/Bagirov_Bartman_Ochal_2023.py +++ b/examples/Bagirov_Bartman_Ochal_2023.py @@ -46,7 +46,7 @@ def potential_normal_direction(u_nu: float) -> float: if u_nu < 1 * mm: return k10 * u_nu ** 2 + k11 * u_nu if u_nu < 2 * mm: - return k20 * u_nu ** 2 + k21 * u_nu + 4.0 + return k20 * u_nu ** 2 + k21 * u_nu + 4 return 16 @staticmethod @@ -95,7 +95,7 @@ def friction_bound(u_nu: float) -> float: def main(save: bool = False): setup = StaticSetup(mesh_type="cross") - for method in ("Powell", "BFGS", "qsm")[2:]: + for method in ("Powell", "BFGS", "CG", "qsm"): for force in np.arange(25e3 * kN, 26e3 * kN + 1, 1e3 * kN) * surface: def outer_forces(x, t=None): if x[1] >= 0.0099: diff --git a/examples/Sofonea_Ochal_Bartman_2023.py b/examples/Sofonea_Ochal_Bartman_2023.py index 306c139d..5f343c54 100644 --- a/examples/Sofonea_Ochal_Bartman_2023.py +++ b/examples/Sofonea_Ochal_Bartman_2023.py @@ -167,6 +167,12 @@ def zero_relaxation(t=None): setup.relaxation = examples[name]["relaxation"] runner = QuasistaticRelaxation(setup, solving_method="schur") + bid = runner.body.mesh.contact_boundary[:, 0] + eid = runner.body.mesh.contact_boundary[:, 1] + bx = runner.body.mesh.initial_nodes[bid][:, 0] + ex = runner.body.mesh.initial_nodes[eid][:, 0] + length = np.max(ex - bx) + print(length) states = runner.solve( n_steps=examples[name]["n_steps"], From d95041ee6125af29563146c2ad6bb1df3f0c739e Mon Sep 17 00:00:00 2001 From: Piotr Date: Wed, 8 Nov 2023 12:57:31 +0100 Subject: [PATCH 06/34] WIP --- examples/Bagirov_Bartman_Ochal_2023.py | 50 ++++++++++++++++---------- 1 file changed, 32 insertions(+), 18 deletions(-) diff --git a/examples/Bagirov_Bartman_Ochal_2023.py b/examples/Bagirov_Bartman_Ochal_2023.py index d3e2ae2e..14f4ba9f 100644 --- a/examples/Bagirov_Bartman_Ochal_2023.py +++ b/examples/Bagirov_Bartman_Ochal_2023.py @@ -7,19 +7,21 @@ from conmech.helpers.config import Config from conmech.mesh.boundaries_description import BoundariesDescription from conmech.plotting.drawer import Drawer -from conmech.scenarios.problems import ContactLaw, Static -from conmech.simulations.problem_solver import Static as StaticProblemSolver +from conmech.properties.mesh_description import RectangleMeshDescription +from conmech.scenarios.problems import ContactLaw, StaticDisplacementProblem +from conmech.simulations.problem_solver import StaticSolver + mesh_density = 4 kN = 1000 mm = 0.001 E = 1.378e8 * kN kappa = 0.3 surface = 5 * mm * 80 * mm -k0 = (30e6 * kN * surface) -k10 = (10e6 * kN * surface) -k11 = (10e3 * kN * surface) -k20 = (5e6 * kN * surface) -k21 = (5e3 * kN * surface) +k0 = 30e6 * kN * surface +k10 = 10e6 * kN * surface +k11 = 10e3 * kN * surface +k20 = 5e6 * kN * surface +k21 = 5e3 * kN * surface def normal_direction(u_nu: float) -> float: @@ -42,11 +44,11 @@ def potential_normal_direction(u_nu: float) -> float: if u_nu <= 0: return 0.0 if u_nu < 0.5 * mm: - return k0 * u_nu ** 2 + return k0 * u_nu**2 if u_nu < 1 * mm: - return k10 * u_nu ** 2 + k11 * u_nu + return k10 * u_nu**2 + k11 * u_nu if u_nu < 2 * mm: - return k20 * u_nu ** 2 + k21 * u_nu + 4 + return k20 * u_nu**2 + k21 * u_nu + 4 return 16 @staticmethod @@ -68,7 +70,7 @@ def regularized_subderivative_tangential_direction( @dataclass() -class StaticSetup(Static): +class StaticSetup(StaticDisplacementProblem): grid_height: ... = 10 * mm elements_number: ... = (mesh_density, 8 * mesh_density) mu_coef: ... = (E * surface) / (2 * (1 + kappa)) @@ -92,11 +94,23 @@ def friction_bound(u_nu: float) -> float: ) -def main(save: bool = False): - setup = StaticSetup(mesh_type="cross") +def main(config: Config): + """ + Entrypoint to example. + + To see result of simulation you need to call from python `main(Config().init())`. + """ + mesh_descr = RectangleMeshDescription( + initial_position=None, + max_element_perimeter=0.25 * 10 * mm, + scale=[2.5 * 10 * mm, 10 * mm], + ) + + setup = StaticSetup(mesh_descr=mesh_descr) for method in ("Powell", "BFGS", "CG", "qsm"): for force in np.arange(25e3 * kN, 26e3 * kN + 1, 1e3 * kN) * surface: + def outer_forces(x, t=None): if x[1] >= 0.0099: return np.array([0, force]) @@ -104,26 +118,26 @@ def outer_forces(x, t=None): setup.outer_forces = outer_forces - runner = StaticProblemSolver(setup, "schur") + runner = StaticSolver(setup, "schur") state = runner.solve( verbose=True, fixed_point_abs_tol=0.001, initial_displacement=setup.initial_displacement, - method=method + method=method, ) - config = Config() drawer = Drawer(state=state, config=config) drawer.colorful = True - drawer.draw(show=not save, save=save, title=f"{method}: {force}") + drawer.draw(show=config.show, save=config.save, title=f"{method}: {force}") if __name__ == "__main__": from matplotlib import pyplot as plt + X = np.linspace(0, -3 * mm, 1000) Y = np.empty(1000) for i in range(1000): Y[i] = MMLV99.potential_normal_direction(X[i]) plt.plot(X, Y) plt.show() - main(save=True) + main(Config().init()) From 74e616b04bb3365613ee71630eb9d15233323cb7 Mon Sep 17 00:00:00 2001 From: Piotr Date: Fri, 10 Nov 2023 10:34:32 +0100 Subject: [PATCH 07/34] draft of subgradient method --- conmech/solvers/optimization/optimization.py | 5 ++--- examples/Bagirov_Bartman_Ochal_2023.py | 8 +++++--- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/conmech/solvers/optimization/optimization.py b/conmech/solvers/optimization/optimization.py index 19f84d24..e4303c7d 100644 --- a/conmech/solvers/optimization/optimization.py +++ b/conmech/solvers/optimization/optimization.py @@ -1,7 +1,6 @@ """ Created at 18.02.2021 """ - import math from typing import Optional @@ -114,9 +113,9 @@ def _solve_impl( "qsmlm", ): # pylint: disable=import-outside-toplevel,import-error) - from kosopt import qsmlm + from kosopt import qsmlmi - solution = qsmlm.minimize(self.loss, solution, args=args, maxiter=maxiter) + solution = qsmlmi.minimize(self.loss, solution, args=args, maxiter=maxiter) sols.append(solution.copy()) elif method.lower() == "constrained": contact_nodes_count = self.body.mesh.boundaries.contact_nodes_count diff --git a/examples/Bagirov_Bartman_Ochal_2023.py b/examples/Bagirov_Bartman_Ochal_2023.py index 14f4ba9f..a09bf002 100644 --- a/examples/Bagirov_Bartman_Ochal_2023.py +++ b/examples/Bagirov_Bartman_Ochal_2023.py @@ -103,13 +103,15 @@ def main(config: Config): mesh_descr = RectangleMeshDescription( initial_position=None, max_element_perimeter=0.25 * 10 * mm, - scale=[2.5 * 10 * mm, 10 * mm], + scale=[8 * 10 * mm, 10 * mm], ) setup = StaticSetup(mesh_descr=mesh_descr) - for method in ("Powell", "BFGS", "CG", "qsm"): - for force in np.arange(25e3 * kN, 26e3 * kN + 1, 1e3 * kN) * surface: + for method in ("Powell", "BFGS", "CG", "qsm")[3:]: + for force in ( + np.asarray([23e3 * kN, 26.2e3 * kN, 27e3 * kN, 30e3 * kN]) * surface + ): def outer_forces(x, t=None): if x[1] >= 0.0099: From bc541998199feb83daba1fc79f88d3b8411a3295 Mon Sep 17 00:00:00 2001 From: Piotr Bartman Date: Mon, 13 Nov 2023 12:36:55 +0100 Subject: [PATCH 08/34] add dg --- conmech/plotting/drawer.py | 3 ++- conmech/solvers/optimization/optimization.py | 10 ++++++++++ conmech/solvers/solver.py | 6 ++++++ examples/Bagirov_Bartman_Ochal_2023.py | 11 ++++++++--- 4 files changed, 26 insertions(+), 4 deletions(-) diff --git a/conmech/plotting/drawer.py b/conmech/plotting/drawer.py index e013325c..6b2f4c41 100644 --- a/conmech/plotting/drawer.py +++ b/conmech/plotting/drawer.py @@ -4,6 +4,7 @@ @author: Michał Jureczka @author: Piotr Bartman """ +import time import matplotlib.pyplot as plt import networkx as nx @@ -252,7 +253,7 @@ def draw_stress(self, axes): def get_output_path(config, format_, name): output_dir = config.output_dir or str(config.timestamp) directory = f"{config.outputs_path}/{output_dir}" - name = name if isinstance(name, str) else config.timestamp + name = name if isinstance(name, str) else time.time_ns() path = f"{directory}/{name}.{format_}" return path diff --git a/conmech/solvers/optimization/optimization.py b/conmech/solvers/optimization/optimization.py index e4303c7d..c90fb12e 100644 --- a/conmech/solvers/optimization/optimization.py +++ b/conmech/solvers/optimization/optimization.py @@ -111,6 +111,16 @@ def _solve_impl( "quasi secant method limited memory", "qsm", "qsmlm", + ): + # pylint: disable=import-outside-toplevel,import-error) + from kosopt import qsmlm + solution = qsmlm.minimize( + self.loss, solution, args=args, maxiter=maxiter + ) + elif method.lower() in ( # TODO + "discontinuous gradient", + "discontinuous gradient method", + "dg", ): # pylint: disable=import-outside-toplevel,import-error) from kosopt import qsmlmi diff --git a/conmech/solvers/solver.py b/conmech/solvers/solver.py index 9b0fdd06..bf4da2a6 100644 --- a/conmech/solvers/solver.py +++ b/conmech/solvers/solver.py @@ -5,6 +5,7 @@ from typing import Optional import numpy as np +import time from conmech.dynamics.statement import Statement, Variables from conmech.dynamics.contact.contact_law import ContactLaw @@ -48,6 +49,8 @@ def __init__( ) ) + self.last_timing = None + def __str__(self) -> str: raise NotImplementedError() @@ -65,6 +68,7 @@ def _solve_impl( raise NotImplementedError() def solve(self, initial_guess: np.ndarray, **kwargs) -> np.ndarray: + start = time.time() solution = self._solve_impl( initial_guess, variable_old=self.v_vector, @@ -80,4 +84,6 @@ def solve(self, initial_guess: np.ndarray, **kwargs) -> np.ndarray: ): solution[i] = c[j] + self.last_timing = time.time() - start + return solution diff --git a/examples/Bagirov_Bartman_Ochal_2023.py b/examples/Bagirov_Bartman_Ochal_2023.py index a09bf002..2db00671 100644 --- a/examples/Bagirov_Bartman_Ochal_2023.py +++ b/examples/Bagirov_Bartman_Ochal_2023.py @@ -108,7 +108,7 @@ def main(config: Config): setup = StaticSetup(mesh_descr=mesh_descr) - for method in ("Powell", "BFGS", "CG", "qsm")[3:]: + for method in ("Powell", "BFGS", "CG", "qsm", "dg")[:]: for force in ( np.asarray([23e3 * kN, 26.2e3 * kN, 27e3 * kN, 30e3 * kN]) * surface ): @@ -130,7 +130,12 @@ def outer_forces(x, t=None): ) drawer = Drawer(state=state, config=config) drawer.colorful = True - drawer.draw(show=config.show, save=config.save, title=f"{method}: {force}") + drawer.draw( + show=config.show, + save=config.save, + title=f"{method}: {force}, " + f"time: {runner.step_solver.last_timing}" + ) if __name__ == "__main__": @@ -142,4 +147,4 @@ def outer_forces(x, t=None): Y[i] = MMLV99.potential_normal_direction(X[i]) plt.plot(X, Y) plt.show() - main(Config().init()) + main(Config(save=True, show=False).init()) From 8a718aeb01b44940ef3a4b2935111a4d5854812a Mon Sep 17 00:00:00 2001 From: Piotr Bartman Date: Sun, 18 Feb 2024 02:20:34 +1100 Subject: [PATCH 09/34] working WIP --- conmech/solvers/optimization/optimization.py | 7 ++ examples/Bagirov_Bartman_Ochal_2023.py | 119 ++++++++++++++++--- 2 files changed, 112 insertions(+), 14 deletions(-) diff --git a/conmech/solvers/optimization/optimization.py b/conmech/solvers/optimization/optimization.py index c90fb12e..fae5b0e4 100644 --- a/conmech/solvers/optimization/optimization.py +++ b/conmech/solvers/optimization/optimization.py @@ -117,6 +117,13 @@ def _solve_impl( solution = qsmlm.minimize( self.loss, solution, args=args, maxiter=maxiter ) + elif method.lower() in ( + "subgradient" + ): + from kosopt import subgradient + solution = subgradient.minimize( + self.loss, solution, args=args, maxiter=maxiter + ) elif method.lower() in ( # TODO "discontinuous gradient", "discontinuous gradient method", diff --git a/examples/Bagirov_Bartman_Ochal_2023.py b/examples/Bagirov_Bartman_Ochal_2023.py index 2db00671..64c42ba3 100644 --- a/examples/Bagirov_Bartman_Ochal_2023.py +++ b/examples/Bagirov_Bartman_Ochal_2023.py @@ -1,6 +1,7 @@ """ Created at 21.08.2019 """ +import pickle from dataclasses import dataclass import numpy as np @@ -10,6 +11,7 @@ from conmech.properties.mesh_description import RectangleMeshDescription from conmech.scenarios.problems import ContactLaw, StaticDisplacementProblem from conmech.simulations.problem_solver import StaticSolver +from conmech.solvers.optimization.optimization import Optimization mesh_density = 4 kN = 1000 @@ -44,11 +46,11 @@ def potential_normal_direction(u_nu: float) -> float: if u_nu <= 0: return 0.0 if u_nu < 0.5 * mm: - return k0 * u_nu**2 + return k0 * u_nu ** 2 if u_nu < 1 * mm: - return k10 * u_nu**2 + k11 * u_nu + return k10 * u_nu ** 2 + k11 * u_nu if u_nu < 2 * mm: - return k20 * u_nu**2 + k21 * u_nu + 4 + return k20 * u_nu ** 2 + k21 * u_nu + 4 return 16 @staticmethod @@ -61,7 +63,7 @@ def subderivative_normal_direction(u_nu: float, v_nu: float) -> float: @staticmethod def regularized_subderivative_tangential_direction( - u_tau: np.ndarray, v_tau: np.ndarray, rho=1e-7 + u_tau: np.ndarray, v_tau: np.ndarray, rho=1e-7 ) -> float: """ Coulomb regularization @@ -94,28 +96,43 @@ def friction_bound(u_nu: float) -> float: ) -def main(config: Config): +def main(config: Config, methods, forces): """ Entrypoint to example. To see result of simulation you need to call from python `main(Config().init())`. """ + PREFIX = "BBO" + if config.force: + to_simulate = [(m, f) for m in methods for f in forces] + else: + to_simulate = [] + for m in methods: + for f in forces: + try: + path = f"{config.outputs_path}/{PREFIX}_mtd_{m}_frc_{f:.2e}" + with open(path, "rb") as output: + _ = pickle.load(output) + except IOError as ioe: + print(ioe) + to_simulate.append((m, f)) + mesh_descr = RectangleMeshDescription( initial_position=None, max_element_perimeter=0.25 * 10 * mm, scale=[8 * 10 * mm, 10 * mm], ) - setup = StaticSetup(mesh_descr=mesh_descr) + if to_simulate: + print("Simulating...") + setup = StaticSetup(mesh_descr=mesh_descr) - for method in ("Powell", "BFGS", "CG", "qsm", "dg")[:]: - for force in ( - np.asarray([23e3 * kN, 26.2e3 * kN, 27e3 * kN, 30e3 * kN]) * surface - ): + for method, force in to_simulate: + print(method, force) def outer_forces(x, t=None): if x[1] >= 0.0099: - return np.array([0, force]) + return np.array([0, force * surface]) return np.array([0, 0]) setup.outer_forces = outer_forces @@ -128,14 +145,38 @@ def outer_forces(x, t=None): initial_displacement=setup.initial_displacement, method=method, ) + path = f"{config.outputs_path}/{PREFIX}_mtd_{method}_frc_{force:.2e}" + with open(path, "wb+") as output: + state.body.dynamics.force.outer.source = None + state.body.dynamics.force.inner.source = None + state.body.properties.relaxation = None + state.setup = None + state.constitutive_law = None + pickle.dump(state, output) + print(Optimization.RESULTS) + + for m in methods: + for f in forces: + path = f"{config.outputs_path}/{PREFIX}_mtd_{m}_frc_{f:.2e}" + with open(path, "rb") as output: + state = pickle.load(output) + drawer = Drawer(state=state, config=config) drawer.colorful = True drawer.draw( show=config.show, save=config.save, - title=f"{method}: {force}, " - f"time: {runner.step_solver.last_timing}" + title=f"{m}: {f}, " + # f"time: {runner.step_solver.last_timing}" ) + x = state.body.mesh.nodes[:state.body.mesh.contact_nodes_count - 1, + 0] + u = state.displacement[:state.body.mesh.contact_nodes_count - 1, 1] + y1 = [normal_direction(-u_) for u_ in u] + plt.plot(x, y1, label=f"{f:.2e}") + plt.title(m) + plt.legend() + plt.show() if __name__ == "__main__": @@ -147,4 +188,54 @@ def outer_forces(x, t=None): Y[i] = MMLV99.potential_normal_direction(X[i]) plt.plot(X, Y) plt.show() - main(Config(save=True, show=False).init()) + # for i in range(1000): + # Y[i] = normal_direction(X[i]) + # plt.plot(X, Y) + # plt.show() + results = { + "Powell": [-0.09786211600599237, + -0.12214289905239312, + -0.13027212877878766, + -0.13447218948364842, + -0.13588717514960513, + -0.1373096435316275, + -0.7582249893948801, + -0.8589012530191608, + -1.2688709210981735, ], + "BFGS": [-0.09487765162385353, + -0.12207326519092926, + -0.11772218280878324, + -0.1198269497911567, + -0.12061095335641955, + -0.1219350781729528, + -0.12279409585312624, + -0.8584230093357013, + -1.2687124265093166, ], + "CG": [-0.0955742828809952, + -0.12191044159984168, + -0.13009806547436803, + -0.1341887930175023, + -0.1358025353476277, + -0.136904521914724, + -0.13865495481609302, + -0.8584104766729636, + -1.2658836730355307, ], + "subgradient2": [-0.09786204500205781, + -0.12214281874382482, + -0.13027204041320914, + -0.15450061948841598, + -0.1571765749815719, + -0.15986547858657962, + -0.7582249071621823, + -0.8589012119331098, + -1.2688708874747643, ], + } + for m, losses in results.items(): + plt.plot(-1 * np.asarray(losses), label=m) + plt.legend() + # plt.loglog() + plt.show() + methods = ("BFGS", "CG", "qsm", "Powell", "subgradient")[-2:] + forces = (23e3 * kN, 25e3 * kN, 25.6e3 * kN, 25.9e3 * kN, 26e3 * kN, + 26.1e3 * kN, 26.2e3 * kN, 27e3 * kN, 30e3 * kN)[:] + main(Config(save=True, show=False, force=False).init(), methods, forces) From 381777f3e0dd3d6358e909eca054b685d556d48b Mon Sep 17 00:00:00 2001 From: Piotr Bartman Date: Tue, 26 Mar 2024 13:52:45 +0100 Subject: [PATCH 10/34] WIP --- examples/Bagirrov_Bartman_Ochal_2024.py | 241 ++++++++++++++++++ ...man_Ochal_2023.py => Makela_et_al_1998.py} | 0 2 files changed, 241 insertions(+) create mode 100644 examples/Bagirrov_Bartman_Ochal_2024.py rename examples/{Bagirov_Bartman_Ochal_2023.py => Makela_et_al_1998.py} (100%) diff --git a/examples/Bagirrov_Bartman_Ochal_2024.py b/examples/Bagirrov_Bartman_Ochal_2024.py new file mode 100644 index 00000000..a880b768 --- /dev/null +++ b/examples/Bagirrov_Bartman_Ochal_2024.py @@ -0,0 +1,241 @@ +""" +Created at 21.08.2019 +""" +import pickle +from dataclasses import dataclass +from matplotlib import pyplot as plt + +import numpy as np +from conmech.helpers.config import Config +from conmech.mesh.boundaries_description import BoundariesDescription +from conmech.plotting.drawer import Drawer +from conmech.properties.mesh_description import RectangleMeshDescription +from conmech.scenarios.problems import ContactLaw, StaticDisplacementProblem +from conmech.simulations.problem_solver import StaticSolver +from conmech.solvers.optimization.optimization import Optimization + +mesh_density = 4 +kN = 1000 +mm = 0.001 +E = 1.378e8 * kN +kappa = 0.3 +surface = 5 * mm * 80 * mm +k0 = 30e6 * kN * surface +k10 = 10e6 * kN * surface +k11 = 10e3 * kN * surface +k20 = 5e6 * kN * surface +k21 = 5e3 * kN * surface +k30 = 2.5e12 * kN * surface + + +def normal_direction(u_nu: float) -> float: + if u_nu <= 0: + return 0.0 + if u_nu < 0.5 * mm: + return k0 * u_nu * 2 + if u_nu < 1 * mm: + return k10 * (u_nu * 2) + k11 + if u_nu < 2 * mm: + return k20 * (u_nu * 2) + k21 + return u_nu ** 3 * 4 * k30 + + +class MMLV99(ContactLaw): + @staticmethod + def potential_normal_direction(u_nu: float) -> float: + if u_nu <= 0: + return 0.0 + if u_nu < 0.5 * mm: + return k0 * u_nu ** 2 + if u_nu < 1 * mm: + return k10 * u_nu ** 2 + k11 * u_nu + if u_nu < 2 * mm: + return k20 * u_nu ** 2 + k21 * u_nu + 4 + return u_nu ** 4 * k30 + + @staticmethod + def potential_tangential_direction(u_tau: np.ndarray) -> float: + return np.log(np.sum(u_tau * u_tau) ** 0.5 + 1) + + @staticmethod + def subderivative_normal_direction(u_nu: float, v_nu: float) -> float: + return 0 + + @staticmethod + def regularized_subderivative_tangential_direction( + u_tau: np.ndarray, v_tau: np.ndarray, rho=1e-7 + ) -> float: + """ + Coulomb regularization + """ + return 0 + + +@dataclass() +class StaticSetup(StaticDisplacementProblem): + grid_height: ... = 10 * mm + elements_number: ... = (mesh_density, 8 * mesh_density) + mu_coef: ... = (E * surface) / (2 * (1 + kappa)) + la_coef: ... = ((E * surface) * kappa) / ((1 + kappa) * (1 - 2 * kappa)) + contact_law: ... = MMLV99 + + @staticmethod + def inner_forces(x, t=None): + return np.array([0.0, 0.0]) + + @staticmethod + def outer_forces(x, t=None): + return np.array([0, 5]) + + @staticmethod + def friction_bound(u_nu: float) -> float: + return 0.0 + + boundaries: ... = BoundariesDescription( + contact=lambda x: x[1] == 0, dirichlet=lambda x: x[0] == 0 + ) + + +def main(config: Config, methods, forces): + """ + Entrypoint to example. + + To see result of simulation you need to call from python `main(Config().init())`. + """ + PREFIX = "DBBO" + if config.force: + to_simulate = [(m, f) for m in methods for f in forces] + else: + to_simulate = [] + for m in methods: + for f in forces: + try: + path = f"{config.outputs_path}/{PREFIX}_mtd_{m}_frc_{f:.2e}" + with open(path, "rb") as output: + _ = pickle.load(output) + except IOError as ioe: + print(ioe) + to_simulate.append((m, f)) + + mesh_descr = RectangleMeshDescription( + initial_position=None, + max_element_perimeter=0.25 * 10 * mm, + scale=[8 * 10 * mm, 10 * mm], + ) + + if to_simulate: + print("Simulating...") + setup = StaticSetup(mesh_descr=mesh_descr) + + for method, force in to_simulate: + print(method, force) + + def outer_forces(x, t=None): + if x[1] >= 0.0099: + return np.array([0, -1 * force * surface]) + return np.array([0, 0]) + + setup.outer_forces = outer_forces + + runner = StaticSolver(setup, "schur") + + state = runner.solve( + verbose=True, + fixed_point_abs_tol=0.001, + initial_displacement=setup.initial_displacement, + method=method, + maxiter=100, + ) + path = f"{config.outputs_path}/{PREFIX}_mtd_{method}_frc_{force:.2e}" + with open(path, "wb+") as output: + state.body.dynamics.force.outer.source = None + state.body.dynamics.force.inner.source = None + state.body.properties.relaxation = None + state.setup = None + state.constitutive_law = None + pickle.dump(state, output) + print(Optimization.RESULTS) + + for m in methods: + for f in forces: + path = f"{config.outputs_path}/{PREFIX}_mtd_{m}_frc_{f:.2e}" + with open(path, "rb") as output: + state = pickle.load(output) + + # drawer = Drawer(state=state, config=config) + # drawer.colorful = True + # drawer.draw( + # show=config.show, + # save=config.save, + # title=f"{m}: {f}, " + # # f"time: {runner.step_solver.last_timing}" + # ) + x = state.body.mesh.nodes[:state.body.mesh.contact_nodes_count - 1, + 0] + u = state.displacement[:state.body.mesh.contact_nodes_count - 1, 1] + y1 = [normal_direction(-u_) for u_ in u] + print(f) + plt.plot(x, y1, label=f"{f:.2e}") + plt.title(m) + plt.legend() + plt.show() + + +if __name__ == "__main__": + X = np.linspace((2-2) * mm, (2+2) * mm, 1000) + Y = np.empty(1000) + for i in range(1000): + Y[i] = MMLV99.potential_normal_direction(X[i]) + plt.plot(X, Y) + plt.show() + for i in range(1000): + Y[i] = normal_direction(X[i]) + plt.plot(X, Y) + plt.show() + # results = { + # "Powell": [-0.09786211600599237, + # -0.12214289905239312, + # -0.13027212877878766, + # -0.13447218948364842, + # -0.13588717514960513, + # -0.1373096435316275, + # -0.7582249893948801, + # -0.8589012530191608, + # -1.2688709210981735, ], + # "BFGS": [-0.09487765162385353, + # -0.12207326519092926, + # -0.11772218280878324, + # -0.1198269497911567, + # -0.12061095335641955, + # -0.1219350781729528, + # -0.12279409585312624, + # -0.8584230093357013, + # -1.2687124265093166, ], + # "CG": [-0.0955742828809952, + # -0.12191044159984168, + # -0.13009806547436803, + # -0.1341887930175023, + # -0.1358025353476277, + # -0.136904521914724, + # -0.13865495481609302, + # -0.8584104766729636, + # -1.2658836730355307, ], + # "subgradient2": [-0.09786204500205781, + # -0.12214281874382482, + # -0.13027204041320914, + # -0.15450061948841598, + # -0.1571765749815719, + # -0.15986547858657962, + # -0.7582249071621823, + # -0.8589012119331098, + # -1.2688708874747643, ], + # } + # for m, losses in results.items(): + # plt.plot(-1 * np.asarray(losses), label=m) + # plt.legend() + # # plt.loglog() + # plt.show() + methods = ("BFGS", "CG", "qsm", "Powell", "subgradient", 'qsm')[-1:] + forces = (23e3 * kN, 25e3 * kN, 25.6e3 * kN, 25.9e3 * kN, 26e3 * kN, + 26.1e3 * kN, 26.2e3 * kN, 27e3 * kN, 30e3 * kN)[-1::4] + main(Config(save=False, show=True, force=True).init(), methods, forces) diff --git a/examples/Bagirov_Bartman_Ochal_2023.py b/examples/Makela_et_al_1998.py similarity index 100% rename from examples/Bagirov_Bartman_Ochal_2023.py rename to examples/Makela_et_al_1998.py From eaacde6bce43a27eb4be098af56ceb3ab6915dc5 Mon Sep 17 00:00:00 2001 From: Piotr Bartman Date: Mon, 17 Jun 2024 10:09:30 +0200 Subject: [PATCH 11/34] BBO: Coruna --- conmech/solvers/optimization/optimization.py | 14 +- conmech/solvers/solver_methods.py | 43 ++++++ examples/Bagirrov_Bartman_Ochal_2024.py | 20 +-- examples/Makela_et_al_1998.py | 131 ++++++++++--------- 4 files changed, 133 insertions(+), 75 deletions(-) diff --git a/conmech/solvers/optimization/optimization.py b/conmech/solvers/optimization/optimization.py index fae5b0e4..b97b2f8e 100644 --- a/conmech/solvers/optimization/optimization.py +++ b/conmech/solvers/optimization/optimization.py @@ -16,7 +16,8 @@ from conmech.dynamics.contact.contact_law import PotentialOfContactLaw from conmech.dynamics.contact.interior_contact_law import InteriorContactLaw from conmech.solvers.solver import Solver -from conmech.solvers.solver_methods import make_cost_functional, make_equation +from conmech.solvers.solver_methods import ( + make_cost_functional, make_equation, make_cost_functional_subgradient) class Optimization(Solver): @@ -43,6 +44,11 @@ def __init__( variable_dimension=statement.dimension_out, problem_dimension=statement.dimension_in, ) + self.subgradient = make_cost_functional_subgradient( + djn=contact_law.subderivative_normal_direction, # TODO + djt=None, + dh_functional=None, + ) if isinstance(statement, WaveStatement): if isinstance(contact_law, InteriorContactLaw): self.loss = make_equation( # TODO! @@ -115,14 +121,16 @@ def _solve_impl( # pylint: disable=import-outside-toplevel,import-error) from kosopt import qsmlm solution = qsmlm.minimize( - self.loss, solution, args=args, maxiter=maxiter + self.loss, solution, args=args, maxiter=maxiter, + subgradient=self.subgradient ) elif method.lower() in ( "subgradient" ): from kosopt import subgradient solution = subgradient.minimize( - self.loss, solution, args=args, maxiter=maxiter + self.loss, solution, args=args, maxiter=maxiter, + subgradient=self.subgradient ) elif method.lower() in ( # TODO "discontinuous gradient", diff --git a/conmech/solvers/solver_methods.py b/conmech/solvers/solver_methods.py index ee763124..73cb9335 100644 --- a/conmech/solvers/solver_methods.py +++ b/conmech/solvers/solver_methods.py @@ -210,3 +210,46 @@ def cost_functional( return result return cost_functional + +def make_cost_functional_subgradient( + djn: Callable, djt: Optional[Callable] = None, dh_functional: Optional[Callable] = None +): + djn = njit(djn) + djt = njit(djt) + dh_functional = njit(dh_functional) + + @numba.njit() + def contact_subgradient(u_vector, u_vector_old, nodes, contact_boundary, + contact_normals): + cost = np.zeros_like(u_vector) + offset = len(u_vector) // DIMENSION + + for edge in contact_boundary: + n_id_0 = edge[0] + n_id_1 = edge[1] + n_0 = nodes[n_id_0] + n_1 = nodes[n_id_1] + if n_id_0 < offset: + um_normal_0 = -n_0[0] # TODO + cost[n_id_0] = djn(um_normal_0) + cost[n_id_0 + offset] = cost[n_id_0] + if n_id_1 < offset: + um_normal_1 = -n_1[0] # TODO + cost[n_id_1] = djn(um_normal_1) + cost[n_id_1 + offset] = cost[n_id_1] + return cost + + # pylint: disable=unused-argument # 'dt' + @numba.njit() + def subgradient( + u_vector, nodes, contact_boundary, contact_normals, lhs, rhs, + u_vector_old, dt + ): + dj = contact_subgradient( + u_vector, u_vector_old, nodes, contact_boundary, contact_normals + ) + result = np.dot(lhs, u_vector) - rhs + dj + result = result.ravel() + return result + + return subgradient diff --git a/examples/Bagirrov_Bartman_Ochal_2024.py b/examples/Bagirrov_Bartman_Ochal_2024.py index a880b768..604eb279 100644 --- a/examples/Bagirrov_Bartman_Ochal_2024.py +++ b/examples/Bagirrov_Bartman_Ochal_2024.py @@ -162,14 +162,14 @@ def outer_forces(x, t=None): with open(path, "rb") as output: state = pickle.load(output) - # drawer = Drawer(state=state, config=config) - # drawer.colorful = True - # drawer.draw( - # show=config.show, - # save=config.save, - # title=f"{m}: {f}, " - # # f"time: {runner.step_solver.last_timing}" - # ) + drawer = Drawer(state=state, config=config) + drawer.colorful = True + drawer.draw( + show=config.show, + save=config.save, + title=f"{m}: {f}, " + # f"time: {runner.step_solver.last_timing}" + ) x = state.body.mesh.nodes[:state.body.mesh.contact_nodes_count - 1, 0] u = state.displacement[:state.body.mesh.contact_nodes_count - 1, 1] @@ -235,7 +235,7 @@ def outer_forces(x, t=None): # plt.legend() # # plt.loglog() # plt.show() - methods = ("BFGS", "CG", "qsm", "Powell", "subgradient", 'qsm')[-1:] + methods = ("BFGS", "CG", "qsm", "Powell", "subgradient") forces = (23e3 * kN, 25e3 * kN, 25.6e3 * kN, 25.9e3 * kN, 26e3 * kN, 26.1e3 * kN, 26.2e3 * kN, 27e3 * kN, 30e3 * kN)[-1::4] - main(Config(save=False, show=True, force=True).init(), methods, forces) + main(Config(save=False, show=True, force=False).init(), methods, forces) diff --git a/examples/Makela_et_al_1998.py b/examples/Makela_et_al_1998.py index 64c42ba3..47525116 100644 --- a/examples/Makela_et_al_1998.py +++ b/examples/Makela_et_al_1998.py @@ -26,32 +26,33 @@ k21 = 5e3 * kN * surface -def normal_direction(u_nu: float) -> float: - u_nu = -u_nu - if u_nu <= 0: - return 0.0 - if u_nu < 0.5 * mm: - return k0 * u_nu * 2 - if u_nu < 1 * mm: - return k10 * (u_nu * 2) + k11 - if u_nu < 2 * mm: - return k20 * (u_nu * 2) + k21 - return 0 - - class MMLV99(ContactLaw): @staticmethod def potential_normal_direction(u_nu: float) -> float: u_nu = -u_nu + coef = 1. if u_nu <= 0: return 0.0 if u_nu < 0.5 * mm: - return k0 * u_nu ** 2 + return k0 * u_nu ** 2 * coef if u_nu < 1 * mm: - return k10 * u_nu ** 2 + k11 * u_nu + return (k10 * u_nu ** 2 + k11 * u_nu) * coef if u_nu < 2 * mm: - return k20 * u_nu ** 2 + k21 * u_nu + 4 - return 16 + return (k20 * u_nu ** 2 + k21 * u_nu + 4) * coef + return 16 * coef + + @staticmethod + def normal_direction(u_nu: float) -> float: + u_nu = -u_nu + if u_nu <= 0: + return 0.0 + if u_nu < 0.5 * mm: + return k0 * u_nu * 2 + if u_nu < 1 * mm: + return k10 * (u_nu * 2) + k11 + if u_nu < 2 * mm: + return k20 * (u_nu * 2) + k21 + return 0 @staticmethod def potential_tangential_direction(u_tau: np.ndarray) -> float: @@ -161,19 +162,23 @@ def outer_forces(x, t=None): with open(path, "rb") as output: state = pickle.load(output) + print("drawing") drawer = Drawer(state=state, config=config) drawer.colorful = True drawer.draw( show=config.show, save=config.save, - title=f"{m}: {f}, " + # title=f"{m}: {f}, " # f"time: {runner.step_solver.last_timing}" ) x = state.body.mesh.nodes[:state.body.mesh.contact_nodes_count - 1, 0] u = state.displacement[:state.body.mesh.contact_nodes_count - 1, 1] - y1 = [normal_direction(-u_) for u_ in u] + y1 = [MMLV99().normal_direction(-u_) for u_ in u] plt.plot(x, y1, label=f"{f:.2e}") + plt.ylabel("Interlaminar binding force [kN/m$^2$]") + plt.xlabel(r"Contact interface [mm]") + plt.grid() plt.title(m) plt.legend() plt.show() @@ -188,54 +193,56 @@ def outer_forces(x, t=None): Y[i] = MMLV99.potential_normal_direction(X[i]) plt.plot(X, Y) plt.show() - # for i in range(1000): - # Y[i] = normal_direction(X[i]) - # plt.plot(X, Y) - # plt.show() + for i in range(1000): + Y[i] = MMLV99.normal_direction(X[i]) + plt.plot(X, Y) + plt.show() results = { - "Powell": [-0.09786211600599237, - -0.12214289905239312, - -0.13027212877878766, - -0.13447218948364842, - -0.13588717514960513, - -0.1373096435316275, - -0.7582249893948801, - -0.8589012530191608, - -1.2688709210981735, ], - "BFGS": [-0.09487765162385353, - -0.12207326519092926, - -0.11772218280878324, - -0.1198269497911567, - -0.12061095335641955, - -0.1219350781729528, - -0.12279409585312624, - -0.8584230093357013, - -1.2687124265093166, ], - "CG": [-0.0955742828809952, + "BFGS": [-0.061546678021737036, + -0.06782602334922566, + -0.07441012406759984, + -0.08129875924227234, + -0.0959892642846613, + -0.10379118250601398, + -0.10538811134540409, + -0.8584224789292736, + -0.14133884664811114, ], + "CG": [-0.07225702623584927, + -0.07966800277816762, + -0.08744039267159345, + -0.09557428287965247, -0.12191044159984168, - -0.13009806547436803, - -0.1341887930175023, -0.1358025353476277, - -0.136904521914724, -0.13865495481609302, - -0.8584104766729636, - -1.2658836730355307, ], - "subgradient2": [-0.09786204500205781, - -0.12214281874382482, - -0.13027204041320914, - -0.15450061948841598, - -0.1571765749815719, - -0.15986547858657962, - -0.7582249071621823, - -0.8589012119331098, - -1.2688708874747643, ], + -0.15028696247286885, + -1.265832916470563, ], + "Powell": [-0.0723012449592487, + -0.07971212256709342, + -0.0874845064006726, + -0.0978621160055679, + -0.12214289905071576, + -0.13588717513833654, + -0.7582249892835198, + -0.8589012526317955, + -1.2688709207679356, ], + "subgradient": [-0.05079652409797247, + -0.046161334145372934, + -0.04120648554585715, + -0.3859157295854724, + -0.6104716467978587, + -0.7302821710666211, + -0.7554950402698594, + -0.8555741662642888, + -1.2663638426265278, ], } - for m, losses in results.items(): - plt.plot(-1 * np.asarray(losses), label=m) + forces = np.asarray((20e3 * kN, 21e3 * kN, 21e3 * kN, 23e3 * kN, + 25e3 * kN, 26e3 * kN, 26.2e3 * kN, 27e3 * kN, 30e3 * kN))[::2] + # for m, losses in results.items(): + # plt.plot(forces/1e3, -1 * np.asarray(losses[:]), "-o", label=m) plt.legend() - # plt.loglog() + plt.ylabel("$-\mathcal{L}(u)$") + plt.xlabel(r"Load [kN/m$^2$]") + plt.grid() plt.show() - methods = ("BFGS", "CG", "qsm", "Powell", "subgradient")[-2:] - forces = (23e3 * kN, 25e3 * kN, 25.6e3 * kN, 25.9e3 * kN, 26e3 * kN, - 26.1e3 * kN, 26.2e3 * kN, 27e3 * kN, 30e3 * kN)[:] + methods = ("BFGS", "CG", "Powell", "subgradient")[-1:] main(Config(save=True, show=False, force=False).init(), methods, forces) From 93f14a01749053ac83a194418ce86e6893ac9923 Mon Sep 17 00:00:00 2001 From: Piotr Bartman Date: Fri, 26 Jul 2024 08:23:00 +0200 Subject: [PATCH 12/34] add test for optimization --- examples/Makela_et_al_1998.py | 116 ++++---- .../regression/test_Makela_et_al_1998.py | 263 ++++++++++++++++++ 2 files changed, 321 insertions(+), 58 deletions(-) create mode 100644 tests/test_conmech/regression/test_Makela_et_al_1998.py diff --git a/examples/Makela_et_al_1998.py b/examples/Makela_et_al_1998.py index 47525116..a1b82eb2 100644 --- a/examples/Makela_et_al_1998.py +++ b/examples/Makela_et_al_1998.py @@ -78,7 +78,7 @@ class StaticSetup(StaticDisplacementProblem): elements_number: ... = (mesh_density, 8 * mesh_density) mu_coef: ... = (E * surface) / (2 * (1 + kappa)) la_coef: ... = ((E * surface) * kappa) / ((1 + kappa) * (1 - 2 * kappa)) - contact_law: ... = MMLV99 + contact_law: ... = MMLV99() @staticmethod def inner_forces(x, t=None): @@ -187,62 +187,62 @@ def outer_forces(x, t=None): if __name__ == "__main__": from matplotlib import pyplot as plt - X = np.linspace(0, -3 * mm, 1000) - Y = np.empty(1000) - for i in range(1000): - Y[i] = MMLV99.potential_normal_direction(X[i]) - plt.plot(X, Y) - plt.show() - for i in range(1000): - Y[i] = MMLV99.normal_direction(X[i]) - plt.plot(X, Y) - plt.show() - results = { - "BFGS": [-0.061546678021737036, - -0.06782602334922566, - -0.07441012406759984, - -0.08129875924227234, - -0.0959892642846613, - -0.10379118250601398, - -0.10538811134540409, - -0.8584224789292736, - -0.14133884664811114, ], - "CG": [-0.07225702623584927, - -0.07966800277816762, - -0.08744039267159345, - -0.09557428287965247, - -0.12191044159984168, - -0.1358025353476277, - -0.13865495481609302, - -0.15028696247286885, - -1.265832916470563, ], - "Powell": [-0.0723012449592487, - -0.07971212256709342, - -0.0874845064006726, - -0.0978621160055679, - -0.12214289905071576, - -0.13588717513833654, - -0.7582249892835198, - -0.8589012526317955, - -1.2688709207679356, ], - "subgradient": [-0.05079652409797247, - -0.046161334145372934, - -0.04120648554585715, - -0.3859157295854724, - -0.6104716467978587, - -0.7302821710666211, - -0.7554950402698594, - -0.8555741662642888, - -1.2663638426265278, ], - } + # X = np.linspace(0, -3 * mm, 1000) + # Y = np.empty(1000) + # for i in range(1000): + # Y[i] = MMLV99.potential_normal_direction(X[i]) + # plt.plot(X, Y) + # plt.show() + # for i in range(1000): + # Y[i] = MMLV99.normal_direction(X[i]) + # plt.plot(X, Y) + # plt.show() + # results = { + # "BFGS": [-0.061546678021737036, + # -0.06782602334922566, + # -0.07441012406759984, + # -0.08129875924227234, + # -0.0959892642846613, + # -0.10379118250601398, + # -0.10538811134540409, + # -0.8584224789292736, + # -0.14133884664811114, ], + # "CG": [-0.07225702623584927, + # -0.07966800277816762, + # -0.08744039267159345, + # -0.09557428287965247, + # -0.12191044159984168, + # -0.1358025353476277, + # -0.13865495481609302, + # -0.15028696247286885, + # -1.265832916470563, ], + # "Powell": [-0.0723012449592487, + # -0.07971212256709342, + # -0.0874845064006726, + # -0.0978621160055679, + # -0.12214289905071576, + # -0.13588717513833654, + # -0.7582249892835198, + # -0.8589012526317955, + # -1.2688709207679356, ], + # "subgradient": [-0.05079652409797247, + # -0.046161334145372934, + # -0.04120648554585715, + # -0.3859157295854724, + # -0.6104716467978587, + # -0.7302821710666211, + # -0.7554950402698594, + # -0.8555741662642888, + # -1.2663638426265278, ], + # } forces = np.asarray((20e3 * kN, 21e3 * kN, 21e3 * kN, 23e3 * kN, 25e3 * kN, 26e3 * kN, 26.2e3 * kN, 27e3 * kN, 30e3 * kN))[::2] - # for m, losses in results.items(): - # plt.plot(forces/1e3, -1 * np.asarray(losses[:]), "-o", label=m) - plt.legend() - plt.ylabel("$-\mathcal{L}(u)$") - plt.xlabel(r"Load [kN/m$^2$]") - plt.grid() - plt.show() - methods = ("BFGS", "CG", "Powell", "subgradient")[-1:] - main(Config(save=True, show=False, force=False).init(), methods, forces) + # # for m, losses in results.items(): + # # plt.plot(forces/1e3, -1 * np.asarray(losses[:]), "-o", label=m) + # plt.legend() + # plt.ylabel("$-\mathcal{L}(u)$") + # plt.xlabel(r"Load [kN/m$^2$]") + # plt.grid() + # plt.show() + methods = ("BFGS", "CG", "Powell", "subgradient")[-2:] + main(Config(save=False, show=True, force=True).init(), methods, forces) diff --git a/tests/test_conmech/regression/test_Makela_et_al_1998.py b/tests/test_conmech/regression/test_Makela_et_al_1998.py new file mode 100644 index 00000000..bacfaef2 --- /dev/null +++ b/tests/test_conmech/regression/test_Makela_et_al_1998.py @@ -0,0 +1,263 @@ +# CONMECH @ Jagiellonian University in Kraków +# +# Copyright (C) 2024 Piotr Bartman-Szwarc +# +# This program is free software; you can redistribute it and/or +# modify it under the terms of the GNU General Public License +# as published by the Free Software Foundation; either version 3 +# of the License, or (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, +# USA. +from dataclasses import dataclass + +import numpy as np +import pytest + +from conmech.mesh.boundaries_description import BoundariesDescription +from conmech.scenarios.problems import StaticDisplacementProblem +from conmech.simulations.problem_solver import StaticSolver +from conmech.properties.mesh_description import CrossMeshDescription, RectangleMeshDescription +from conmech.solvers.optimization.optimization import Optimization +from examples.Makela_et_al_1998 import MMLV99 +from tests.test_conmech.regression.std_boundary import standard_boundary_nodes + + +@pytest.fixture(params=["schur"]) +def solving_method(request): + return request.param + + +def generate_test_suits(): + test_suites = [] + + optimization_mtd_pow = "Powell" + + expected_displacement_vector_pow = [ + [0., 0.], + [-0.00001052, -0.00000833], + [-0.00001962, -0.00001745], + [-0.00002745, -0.00003174], + [-0.0000345, -0.0000481], + [-0.00004035, -0.00006865], + [-0.00004552, -0.00009066], + [-0.00004975, -0.0001157], + [-0.00005346, -0.00014177], + [-0.00005645, -0.00017], + [-0.00005907, -0.00019895], + [-0.00006118, -0.00022943], + [-0.00006303, -0.00026042], + [-0.00006451, -0.00029252], + [-0.00006582, -0.000325], + [-0.00006687, -0.00035831], + [-0.00006781, -0.00039192], + [-0.00006855, -0.00042619], + [-0.00006919, -0.0004607], + [-0.00006967, -0.00049582], + [-0.00007003, -0.00053111], + [-0.00007025, -0.00056658], + [-0.00007041, -0.0006021], + [-0.00007047, -0.00063776], + [-0.00007051, -0.00067346], + [-0.00007049, -0.00070925], + [-0.00007046, -0.00074507], + [-0.00007041, -0.00078096], + [-0.00007035, -0.00081689], + [-0.00007028, -0.0008529], + [-0.00007019, -0.00088894], + [-0.00007006, -0.00092505], + [-0.00006991, -0.00096119], + [-0.00006972, -0.00099732], + [-0.00003993, -0.00099771], + [-0.00001013, -0.00099814], + [0.00001963, -0.00099859], + [0.00004954, -0.00099912], + [0.00007947, -0.0009997], + [0.00007922, -0.00096314], + [0.00007891, -0.00092685], + [0.00007862, -0.00089072], + [0.00007832, -0.00085463], + [0.00007802, -0.00081863], + [0.00007773, -0.00078266], + [0.00007742, -0.00074677], + [0.00007711, -0.0007109], + [0.00007677, -0.00067512], + [0.00007641, -0.00063936], + [0.00007599, -0.00060371], + [0.00007552, -0.00056812], + [0.00007493, -0.00053273], + [0.00007425, -0.00049744], + [0.00007341, -0.00046252], + [0.00007246, -0.00042781], + [0.00007135, -0.00039366], + [0.00007013, -0.00035978], + [0.00006871, -0.00032663], + [0.00006715, -0.00029382], + [0.00006531, -0.00026197], + [0.00006325, -0.00023055], + [0.00006079, -0.00020043], + [0.00005801, -0.00017093], + [0.00005465, -0.00014324], + [0.00005085, -0.00011642], + [0.00004624, -0.00009215], + [0.00004105, -0.00006911], + [0.00003479, -0.00004963], + [0.00002775, -0.00003184], + [0.00001932, -0.00001891], + [0.00000961, -0.0000077], + [0., 0.], + [0., 0.], + [0., 0.], + [0., 0.], + [0., 0.], + ] + test_suites.append((optimization_mtd_pow, expected_displacement_vector_pow)) + + optimization_mtd_subg = "subgradient" + expected_displacement_vector_subg = [ + [ 0. , 0. ], + [-0.00006594, -0.00004315], + [-0.00012749, -0.00009338], + [-0.00018649, -0.00018164], + [-0.00024325, -0.0002873 ], + [-0.00029591, -0.00043023], + [-0.00034592, -0.00058795], + [-0.00039183, -0.00077834], + [-0.00043521, -0.00098149], + [-0.00047474, -0.00121284], + [-0.00051191, -0.00145513], + [-0.0005455 , -0.00172143], + [-0.0005769 , -0.00199696], + [-0.00060501, -0.00229263], + [-0.00063112, -0.00259596], + [-0.00065425, -0.00291588], + [-0.00067557, -0.00324204], + [-0.00069422, -0.00358162], + [-0.00071126, -0.00392617], + [-0.00072592, -0.00428129], + [-0.00073916, -0.00464026], + [-0.00075031, -0.0050073 ], + [-0.00076022, -0.0053772 ], + [-0.00076831, -0.00575298], + [-0.00077534, -0.00613077], + [-0.00078084, -0.00651256], + [-0.00078546, -0.00689564], + [-0.00078884, -0.00728117], + [-0.00079152, -0.00766739], + [-0.00079327, -0.00805485], + [-0.00079451, -0.00844256], + [-0.00079513, -0.00883062], + [-0.00079547, -0.00921861], + [-0.0007956 , -0.00960656], + [-0.00047553, -0.00960645], + [-0.00015576, -0.00960649], + [ 0.0001639 , -0.00960666], + [ 0.00048371, -0.00960712], + [ 0.0008037 , -0.00960768], + [ 0.00080342, -0.0092195 ], + [ 0.00080296, -0.00883138], + [ 0.0008021 , -0.00844342], + [ 0.00080088, -0.00805553], + [ 0.00079888, -0.0076684 ], + [ 0.0007963 , -0.00728171], + [ 0.00079264, -0.00689683], + [ 0.00078819, -0.00651292], + [ 0.00078235, -0.0061322 ], + [ 0.00077554, -0.00575311], + [ 0.00076705, -0.00537893], + [ 0.00075739, -0.00500716], + [ 0.00074577, -0.00464232], + [ 0.00073279, -0.00428083], + [ 0.00071756, -0.00392862], + [ 0.00070079, -0.00358081], + [ 0.00068145, -0.00324495], + [ 0.00066038, -0.00291468], + [ 0.00063644, -0.00259938], + [ 0.00061058, -0.00229101], + [ 0.00058154, -0.00200098], + [ 0.00055038, -0.00171937], + [ 0.00051574, -0.0014598 ], + [ 0.00047881, -0.0012103 ], + [ 0.00043813, -0.00098686], + [ 0.00039498, -0.00077525], + [ 0.00034784, -0.00059401], + [ 0.00029808, -0.00042646], + [ 0.00024407, -0.00029379], + [ 0.00018721, -0.00017669], + [ 0.00012562, -0.00009923], + [ 0.00005962, -0.00003655], + [ 0. , 0. ], + [ 0. , 0. ], + [ 0. , 0. ], + [ 0. , 0. ], + [ 0. , 0. ], + ] + + test_suites.append((optimization_mtd_subg, expected_displacement_vector_subg)) + + return test_suites + + +@pytest.mark.parametrize("optimization_method, expected_displacement_vector", generate_test_suits()) +def test_static_solver(solving_method, optimization_method, expected_displacement_vector): + mesh_density = 4 + kN = 1000 + mm = 0.001 + E = 1.378e8 * kN + kappa = 0.3 + surface = 5 * mm * 80 * mm + + @dataclass() + class StaticSetup(StaticDisplacementProblem): + grid_height: ... = 10 * mm + elements_number: ... = (mesh_density, 8 * mesh_density) + mu_coef: ... = (E * surface) / (2 * (1 + kappa)) + la_coef: ... = ((E * surface) * kappa) / ((1 + kappa) * (1 - 2 * kappa)) + contact_law: ... = MMLV99() + + @staticmethod + def inner_forces(x, t=None): + return np.array([0.0, 0.0]) + + @staticmethod + def outer_forces(x, t=None): + if x[1] >= 0.0099: + return np.array([0, 26.2e3 * kN * surface]) + return np.array([0, 0]) + + @staticmethod + def friction_bound(u_nu: float) -> float: + return 0.0 + + boundaries: ... = BoundariesDescription( + contact=lambda x: x[1] == 0, dirichlet=lambda x: x[0] == 0 + ) + + mesh_descr = RectangleMeshDescription( + initial_position=None, + max_element_perimeter=0.25 * 10 * mm, + scale=[8 * 10 * mm, 10 * mm], + ) + + setup = StaticSetup(mesh_descr) + + runner = StaticSolver(setup, solving_method) + result = runner.solve(initial_displacement=setup.initial_displacement, method=optimization_method) + + displacement = result.body.mesh.nodes[:] - result.displaced_nodes[:] + std_ids = standard_boundary_nodes(runner.body.mesh.nodes, runner.body.mesh.elements) + + # print result + np.set_printoptions(precision=8, suppress=True) + print(repr(displacement[std_ids])) + + np.testing.assert_array_almost_equal( + displacement[std_ids], expected_displacement_vector, decimal=3 + ) From 544a6de72b14251a9996cbd786a91606d78a9058 Mon Sep 17 00:00:00 2001 From: Piotr Bartman Date: Fri, 26 Jul 2024 10:06:14 +0200 Subject: [PATCH 13/34] working subgradient --- conmech/solvers/optimization/optimization.py | 11 ++++--- conmech/solvers/solver_methods.py | 21 +++++++++---- examples/Makela_et_al_1998.py | 33 ++++++++------------ examples/Sofonea_Ochal_Bartman_2023.py | 4 +-- 4 files changed, 36 insertions(+), 33 deletions(-) diff --git a/conmech/solvers/optimization/optimization.py b/conmech/solvers/optimization/optimization.py index b97b2f8e..bb7bcb7e 100644 --- a/conmech/solvers/optimization/optimization.py +++ b/conmech/solvers/optimization/optimization.py @@ -44,11 +44,12 @@ def __init__( variable_dimension=statement.dimension_out, problem_dimension=statement.dimension_in, ) - self.subgradient = make_cost_functional_subgradient( - djn=contact_law.subderivative_normal_direction, # TODO - djt=None, - dh_functional=None, - ) + if hasattr(contact_law, "subderivative_normal_direction"): + self.subgradient = make_cost_functional_subgradient( + djn=contact_law.subderivative_normal_direction, # TODO + djt=None, + dh_functional=None, + ) if isinstance(statement, WaveStatement): if isinstance(contact_law, InteriorContactLaw): self.loss = make_equation( # TODO! diff --git a/conmech/solvers/solver_methods.py b/conmech/solvers/solver_methods.py index 73cb9335..fd3e7214 100644 --- a/conmech/solvers/solver_methods.py +++ b/conmech/solvers/solver_methods.py @@ -217,6 +217,7 @@ def make_cost_functional_subgradient( djn = njit(djn) djt = njit(djt) dh_functional = njit(dh_functional) + DIMENSION = 2 @numba.njit() def contact_subgradient(u_vector, u_vector_old, nodes, contact_boundary, @@ -231,24 +232,32 @@ def contact_subgradient(u_vector, u_vector_old, nodes, contact_boundary, n_1 = nodes[n_id_1] if n_id_0 < offset: um_normal_0 = -n_0[0] # TODO - cost[n_id_0] = djn(um_normal_0) + cost[n_id_0] = djn(um_normal_0, 0., 0.) cost[n_id_0 + offset] = cost[n_id_0] if n_id_1 < offset: um_normal_1 = -n_1[0] # TODO - cost[n_id_1] = djn(um_normal_1) + cost[n_id_1] = djn(um_normal_1, 0., 0.) cost[n_id_1 + offset] = cost[n_id_1] return cost # pylint: disable=unused-argument # 'dt' @numba.njit() def subgradient( - u_vector, nodes, contact_boundary, contact_normals, lhs, rhs, - u_vector_old, dt + var, + var_old, + nodes, + contact_boundary, + contact_normals, + lhs, + rhs, + u_vector, + base_integrals, + dt, ): dj = contact_subgradient( - u_vector, u_vector_old, nodes, contact_boundary, contact_normals + var, var_old, nodes, contact_boundary, contact_normals ) - result = np.dot(lhs, u_vector) - rhs + dj + result = np.dot(lhs, var) - rhs + dj result = result.ravel() return result diff --git a/examples/Makela_et_al_1998.py b/examples/Makela_et_al_1998.py index a1b82eb2..8b131ae2 100644 --- a/examples/Makela_et_al_1998.py +++ b/examples/Makela_et_al_1998.py @@ -5,11 +5,13 @@ from dataclasses import dataclass import numpy as np + +from conmech.dynamics.contact.contact_law import PotentialOfContactLaw, DirectContactLaw from conmech.helpers.config import Config from conmech.mesh.boundaries_description import BoundariesDescription from conmech.plotting.drawer import Drawer from conmech.properties.mesh_description import RectangleMeshDescription -from conmech.scenarios.problems import ContactLaw, StaticDisplacementProblem +from conmech.scenarios.problems import StaticDisplacementProblem from conmech.simulations.problem_solver import StaticSolver from conmech.solvers.optimization.optimization import Optimization @@ -26,10 +28,10 @@ k21 = 5e3 * kN * surface -class MMLV99(ContactLaw): +class MMLV99(PotentialOfContactLaw, DirectContactLaw): @staticmethod - def potential_normal_direction(u_nu: float) -> float: - u_nu = -u_nu + def potential_normal_direction(var_nu: float, static_displacement_nu: float, dt: float) -> float: + u_nu = -var_nu coef = 1. if u_nu <= 0: return 0.0 @@ -42,8 +44,10 @@ def potential_normal_direction(u_nu: float) -> float: return 16 * coef @staticmethod - def normal_direction(u_nu: float) -> float: - u_nu = -u_nu + def subderivative_normal_direction( + var_nu: float, static_displacement_nu: float, dt: float + ) -> float: + u_nu = -var_nu if u_nu <= 0: return 0.0 if u_nu < 0.5 * mm: @@ -55,21 +59,10 @@ def normal_direction(u_nu: float) -> float: return 0 @staticmethod - def potential_tangential_direction(u_tau: np.ndarray) -> float: - return np.log(np.sum(u_tau * u_tau) ** 0.5 + 1) - - @staticmethod - def subderivative_normal_direction(u_nu: float, v_nu: float) -> float: - return 0 - - @staticmethod - def regularized_subderivative_tangential_direction( - u_tau: np.ndarray, v_tau: np.ndarray, rho=1e-7 + def potential_tangential_direction( + var_tau: float, static_displacement_tau: float, dt: float ) -> float: - """ - Coulomb regularization - """ - return 0 + return np.log(np.sum(var_tau ** 2) ** 0.5 + 1) @dataclass() diff --git a/examples/Sofonea_Ochal_Bartman_2023.py b/examples/Sofonea_Ochal_Bartman_2023.py index 5f343c54..4d83a53a 100644 --- a/examples/Sofonea_Ochal_Bartman_2023.py +++ b/examples/Sofonea_Ochal_Bartman_2023.py @@ -169,8 +169,8 @@ def zero_relaxation(t=None): runner = QuasistaticRelaxation(setup, solving_method="schur") bid = runner.body.mesh.contact_boundary[:, 0] eid = runner.body.mesh.contact_boundary[:, 1] - bx = runner.body.mesh.initial_nodes[bid][:, 0] - ex = runner.body.mesh.initial_nodes[eid][:, 0] + bx = runner.body.mesh.nodes[bid][:, 0] + ex = runner.body.mesh.nodes[eid][:, 0] length = np.max(ex - bx) print(length) From f030e6588b40c8db144e87dbf2fe9985fd0b54eb Mon Sep 17 00:00:00 2001 From: Piotr Bartman-Szwarc Date: Fri, 26 Jul 2024 10:31:33 +0200 Subject: [PATCH 14/34] cleanup --- conmech/plotting/drawer.py | 1 + conmech/solvers/optimization/optimization.py | 45 +++-- conmech/solvers/solver.py | 24 ++- conmech/solvers/solver_methods.py | 34 ++-- examples/Bagirrov_Bartman_Ochal_2024.py | 35 ++-- examples/Makela_et_al_1998.py | 55 ++++-- .../regression/test_Makela_et_al_1998.py | 168 +++++++++--------- 7 files changed, 222 insertions(+), 140 deletions(-) diff --git a/conmech/plotting/drawer.py b/conmech/plotting/drawer.py index 6b2f4c41..43af9857 100644 --- a/conmech/plotting/drawer.py +++ b/conmech/plotting/drawer.py @@ -4,6 +4,7 @@ @author: Michał Jureczka @author: Piotr Bartman """ + import time import matplotlib.pyplot as plt diff --git a/conmech/solvers/optimization/optimization.py b/conmech/solvers/optimization/optimization.py index bb7bcb7e..76e4e857 100644 --- a/conmech/solvers/optimization/optimization.py +++ b/conmech/solvers/optimization/optimization.py @@ -1,6 +1,21 @@ -""" -Created at 18.02.2021 -""" +# CONMECH @ Jagiellonian University in Kraków +# +# Copyright (C) 2021-2024 Piotr Bartman-Szwarc +# +# This program is free software; you can redistribute it and/or +# modify it under the terms of the GNU General Public License +# as published by the Free Software Foundation; either version 3 +# of the License, or (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, +# USA. import math from typing import Optional @@ -17,7 +32,10 @@ from conmech.dynamics.contact.interior_contact_law import InteriorContactLaw from conmech.solvers.solver import Solver from conmech.solvers.solver_methods import ( - make_cost_functional, make_equation, make_cost_functional_subgradient) + make_cost_functional, + make_equation, + make_cost_functional_subgradient, +) class Optimization(Solver): @@ -121,22 +139,21 @@ def _solve_impl( ): # pylint: disable=import-outside-toplevel,import-error) from kosopt import qsmlm + solution = qsmlm.minimize( - self.loss, solution, args=args, maxiter=maxiter, - subgradient=self.subgradient + self.loss, solution, args=args, maxiter=maxiter, subgradient=self.subgradient ) - elif method.lower() in ( - "subgradient" - ): + elif method.lower() in ("subgradient",): + # pylint: disable=import-outside-toplevel,import-error) from kosopt import subgradient + solution = subgradient.minimize( - self.loss, solution, args=args, maxiter=maxiter, - subgradient=self.subgradient + self.loss, solution, args=args, maxiter=maxiter, subgradient=self.subgradient ) elif method.lower() in ( # TODO - "discontinuous gradient", - "discontinuous gradient method", - "dg", + "discontinuous gradient", + "discontinuous gradient method", + "dg", ): # pylint: disable=import-outside-toplevel,import-error) from kosopt import qsmlmi diff --git a/conmech/solvers/solver.py b/conmech/solvers/solver.py index bf4da2a6..296278c4 100644 --- a/conmech/solvers/solver.py +++ b/conmech/solvers/solver.py @@ -1,11 +1,25 @@ -""" -Created at 18.02.2021 -""" - +# CONMECH @ Jagiellonian University in Kraków +# +# Copyright (C) 2021-2024 Piotr Bartman-Szwarc +# +# This program is free software; you can redistribute it and/or +# modify it under the terms of the GNU General Public License +# as published by the Free Software Foundation; either version 3 +# of the License, or (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, +# USA. +import time from typing import Optional import numpy as np -import time from conmech.dynamics.statement import Statement, Variables from conmech.dynamics.contact.contact_law import ContactLaw diff --git a/conmech/solvers/solver_methods.py b/conmech/solvers/solver_methods.py index fd3e7214..f28f15f5 100644 --- a/conmech/solvers/solver_methods.py +++ b/conmech/solvers/solver_methods.py @@ -1,7 +1,21 @@ -""" -Created at 21.08.2019 -""" - +# CONMECH @ Jagiellonian University in Kraków +# +# Copyright (C) 2019-2024 Piotr Bartman-Szwarc +# +# This program is free software; you can redistribute it and/or +# modify it under the terms of the GNU General Public License +# as published by the Free Software Foundation; either version 3 +# of the License, or (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, +# USA. from typing import Callable, Optional, Any import numba @@ -211,6 +225,7 @@ def cost_functional( return cost_functional + def make_cost_functional_subgradient( djn: Callable, djt: Optional[Callable] = None, dh_functional: Optional[Callable] = None ): @@ -220,8 +235,7 @@ def make_cost_functional_subgradient( DIMENSION = 2 @numba.njit() - def contact_subgradient(u_vector, u_vector_old, nodes, contact_boundary, - contact_normals): + def contact_subgradient(u_vector, u_vector_old, nodes, contact_boundary, contact_normals): cost = np.zeros_like(u_vector) offset = len(u_vector) // DIMENSION @@ -232,11 +246,11 @@ def contact_subgradient(u_vector, u_vector_old, nodes, contact_boundary, n_1 = nodes[n_id_1] if n_id_0 < offset: um_normal_0 = -n_0[0] # TODO - cost[n_id_0] = djn(um_normal_0, 0., 0.) + cost[n_id_0] = djn(um_normal_0, 0.0, 0.0) cost[n_id_0 + offset] = cost[n_id_0] if n_id_1 < offset: um_normal_1 = -n_1[0] # TODO - cost[n_id_1] = djn(um_normal_1, 0., 0.) + cost[n_id_1] = djn(um_normal_1, 0.0, 0.0) cost[n_id_1 + offset] = cost[n_id_1] return cost @@ -254,9 +268,7 @@ def subgradient( base_integrals, dt, ): - dj = contact_subgradient( - var, var_old, nodes, contact_boundary, contact_normals - ) + dj = contact_subgradient(var, var_old, nodes, contact_boundary, contact_normals) result = np.dot(lhs, var) - rhs + dj result = result.ravel() return result diff --git a/examples/Bagirrov_Bartman_Ochal_2024.py b/examples/Bagirrov_Bartman_Ochal_2024.py index 604eb279..916dbafd 100644 --- a/examples/Bagirrov_Bartman_Ochal_2024.py +++ b/examples/Bagirrov_Bartman_Ochal_2024.py @@ -1,6 +1,7 @@ """ Created at 21.08.2019 """ + import pickle from dataclasses import dataclass from matplotlib import pyplot as plt @@ -37,7 +38,7 @@ def normal_direction(u_nu: float) -> float: return k10 * (u_nu * 2) + k11 if u_nu < 2 * mm: return k20 * (u_nu * 2) + k21 - return u_nu ** 3 * 4 * k30 + return u_nu**3 * 4 * k30 class MMLV99(ContactLaw): @@ -46,12 +47,12 @@ def potential_normal_direction(u_nu: float) -> float: if u_nu <= 0: return 0.0 if u_nu < 0.5 * mm: - return k0 * u_nu ** 2 + return k0 * u_nu**2 if u_nu < 1 * mm: - return k10 * u_nu ** 2 + k11 * u_nu + return k10 * u_nu**2 + k11 * u_nu if u_nu < 2 * mm: - return k20 * u_nu ** 2 + k21 * u_nu + 4 - return u_nu ** 4 * k30 + return k20 * u_nu**2 + k21 * u_nu + 4 + return u_nu**4 * k30 @staticmethod def potential_tangential_direction(u_tau: np.ndarray) -> float: @@ -63,7 +64,7 @@ def subderivative_normal_direction(u_nu: float, v_nu: float) -> float: @staticmethod def regularized_subderivative_tangential_direction( - u_tau: np.ndarray, v_tau: np.ndarray, rho=1e-7 + u_tau: np.ndarray, v_tau: np.ndarray, rho=1e-7 ) -> float: """ Coulomb regularization @@ -167,12 +168,11 @@ def outer_forces(x, t=None): drawer.draw( show=config.show, save=config.save, - title=f"{m}: {f}, " + title=f"{m}: {f}, ", # f"time: {runner.step_solver.last_timing}" ) - x = state.body.mesh.nodes[:state.body.mesh.contact_nodes_count - 1, - 0] - u = state.displacement[:state.body.mesh.contact_nodes_count - 1, 1] + x = state.body.mesh.nodes[: state.body.mesh.contact_nodes_count - 1, 0] + u = state.displacement[: state.body.mesh.contact_nodes_count - 1, 1] y1 = [normal_direction(-u_) for u_ in u] print(f) plt.plot(x, y1, label=f"{f:.2e}") @@ -182,7 +182,7 @@ def outer_forces(x, t=None): if __name__ == "__main__": - X = np.linspace((2-2) * mm, (2+2) * mm, 1000) + X = np.linspace((2 - 2) * mm, (2 + 2) * mm, 1000) Y = np.empty(1000) for i in range(1000): Y[i] = MMLV99.potential_normal_direction(X[i]) @@ -236,6 +236,15 @@ def outer_forces(x, t=None): # # plt.loglog() # plt.show() methods = ("BFGS", "CG", "qsm", "Powell", "subgradient") - forces = (23e3 * kN, 25e3 * kN, 25.6e3 * kN, 25.9e3 * kN, 26e3 * kN, - 26.1e3 * kN, 26.2e3 * kN, 27e3 * kN, 30e3 * kN)[-1::4] + forces = ( + 23e3 * kN, + 25e3 * kN, + 25.6e3 * kN, + 25.9e3 * kN, + 26e3 * kN, + 26.1e3 * kN, + 26.2e3 * kN, + 27e3 * kN, + 30e3 * kN, + )[-1::4] main(Config(save=False, show=True, force=False).init(), methods, forces) diff --git a/examples/Makela_et_al_1998.py b/examples/Makela_et_al_1998.py index 8b131ae2..070f0998 100644 --- a/examples/Makela_et_al_1998.py +++ b/examples/Makela_et_al_1998.py @@ -1,6 +1,21 @@ -""" -Created at 21.08.2019 -""" +# CONMECH @ Jagiellonian University in Kraków +# +# Copyright (C) 2023-2024 Piotr Bartman-Szwarc +# +# This program is free software; you can redistribute it and/or +# modify it under the terms of the GNU General Public License +# as published by the Free Software Foundation; either version 3 +# of the License, or (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, +# USA. import pickle from dataclasses import dataclass @@ -30,17 +45,19 @@ class MMLV99(PotentialOfContactLaw, DirectContactLaw): @staticmethod - def potential_normal_direction(var_nu: float, static_displacement_nu: float, dt: float) -> float: + def potential_normal_direction( + var_nu: float, static_displacement_nu: float, dt: float + ) -> float: u_nu = -var_nu - coef = 1. + coef = 1.0 if u_nu <= 0: return 0.0 if u_nu < 0.5 * mm: - return k0 * u_nu ** 2 * coef + return k0 * u_nu**2 * coef if u_nu < 1 * mm: - return (k10 * u_nu ** 2 + k11 * u_nu) * coef + return (k10 * u_nu**2 + k11 * u_nu) * coef if u_nu < 2 * mm: - return (k20 * u_nu ** 2 + k21 * u_nu + 4) * coef + return (k20 * u_nu**2 + k21 * u_nu + 4) * coef return 16 * coef @staticmethod @@ -62,7 +79,7 @@ def subderivative_normal_direction( def potential_tangential_direction( var_tau: float, static_displacement_tau: float, dt: float ) -> float: - return np.log(np.sum(var_tau ** 2) ** 0.5 + 1) + return np.log(np.sum(var_tau**2) ** 0.5 + 1) @dataclass() @@ -164,9 +181,8 @@ def outer_forces(x, t=None): # title=f"{m}: {f}, " # f"time: {runner.step_solver.last_timing}" ) - x = state.body.mesh.nodes[:state.body.mesh.contact_nodes_count - 1, - 0] - u = state.displacement[:state.body.mesh.contact_nodes_count - 1, 1] + x = state.body.mesh.nodes[: state.body.mesh.contact_nodes_count - 1, 0] + u = state.displacement[: state.body.mesh.contact_nodes_count - 1, 1] y1 = [MMLV99().normal_direction(-u_) for u_ in u] plt.plot(x, y1, label=f"{f:.2e}") plt.ylabel("Interlaminar binding force [kN/m$^2$]") @@ -228,8 +244,19 @@ def outer_forces(x, t=None): # -0.8555741662642888, # -1.2663638426265278, ], # } - forces = np.asarray((20e3 * kN, 21e3 * kN, 21e3 * kN, 23e3 * kN, - 25e3 * kN, 26e3 * kN, 26.2e3 * kN, 27e3 * kN, 30e3 * kN))[::2] + forces = np.asarray( + ( + 20e3 * kN, + 21e3 * kN, + 21e3 * kN, + 23e3 * kN, + 25e3 * kN, + 26e3 * kN, + 26.2e3 * kN, + 27e3 * kN, + 30e3 * kN, + ) + )[::2] # # for m, losses in results.items(): # # plt.plot(forces/1e3, -1 * np.asarray(losses[:]), "-o", label=m) # plt.legend() diff --git a/tests/test_conmech/regression/test_Makela_et_al_1998.py b/tests/test_conmech/regression/test_Makela_et_al_1998.py index bacfaef2..e6d29e1c 100644 --- a/tests/test_conmech/regression/test_Makela_et_al_1998.py +++ b/tests/test_conmech/regression/test_Makela_et_al_1998.py @@ -41,7 +41,7 @@ def generate_test_suits(): optimization_mtd_pow = "Powell" expected_displacement_vector_pow = [ - [0., 0.], + [0.0, 0.0], [-0.00001052, -0.00000833], [-0.00001962, -0.00001745], [-0.00002745, -0.00003174], @@ -112,92 +112,92 @@ def generate_test_suits(): [0.00002775, -0.00003184], [0.00001932, -0.00001891], [0.00000961, -0.0000077], - [0., 0.], - [0., 0.], - [0., 0.], - [0., 0.], - [0., 0.], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], ] test_suites.append((optimization_mtd_pow, expected_displacement_vector_pow)) optimization_mtd_subg = "subgradient" expected_displacement_vector_subg = [ - [ 0. , 0. ], - [-0.00006594, -0.00004315], - [-0.00012749, -0.00009338], - [-0.00018649, -0.00018164], - [-0.00024325, -0.0002873 ], - [-0.00029591, -0.00043023], - [-0.00034592, -0.00058795], - [-0.00039183, -0.00077834], - [-0.00043521, -0.00098149], - [-0.00047474, -0.00121284], - [-0.00051191, -0.00145513], - [-0.0005455 , -0.00172143], - [-0.0005769 , -0.00199696], - [-0.00060501, -0.00229263], - [-0.00063112, -0.00259596], - [-0.00065425, -0.00291588], - [-0.00067557, -0.00324204], - [-0.00069422, -0.00358162], - [-0.00071126, -0.00392617], - [-0.00072592, -0.00428129], - [-0.00073916, -0.00464026], - [-0.00075031, -0.0050073 ], - [-0.00076022, -0.0053772 ], - [-0.00076831, -0.00575298], - [-0.00077534, -0.00613077], - [-0.00078084, -0.00651256], - [-0.00078546, -0.00689564], - [-0.00078884, -0.00728117], - [-0.00079152, -0.00766739], - [-0.00079327, -0.00805485], - [-0.00079451, -0.00844256], - [-0.00079513, -0.00883062], - [-0.00079547, -0.00921861], - [-0.0007956 , -0.00960656], - [-0.00047553, -0.00960645], - [-0.00015576, -0.00960649], - [ 0.0001639 , -0.00960666], - [ 0.00048371, -0.00960712], - [ 0.0008037 , -0.00960768], - [ 0.00080342, -0.0092195 ], - [ 0.00080296, -0.00883138], - [ 0.0008021 , -0.00844342], - [ 0.00080088, -0.00805553], - [ 0.00079888, -0.0076684 ], - [ 0.0007963 , -0.00728171], - [ 0.00079264, -0.00689683], - [ 0.00078819, -0.00651292], - [ 0.00078235, -0.0061322 ], - [ 0.00077554, -0.00575311], - [ 0.00076705, -0.00537893], - [ 0.00075739, -0.00500716], - [ 0.00074577, -0.00464232], - [ 0.00073279, -0.00428083], - [ 0.00071756, -0.00392862], - [ 0.00070079, -0.00358081], - [ 0.00068145, -0.00324495], - [ 0.00066038, -0.00291468], - [ 0.00063644, -0.00259938], - [ 0.00061058, -0.00229101], - [ 0.00058154, -0.00200098], - [ 0.00055038, -0.00171937], - [ 0.00051574, -0.0014598 ], - [ 0.00047881, -0.0012103 ], - [ 0.00043813, -0.00098686], - [ 0.00039498, -0.00077525], - [ 0.00034784, -0.00059401], - [ 0.00029808, -0.00042646], - [ 0.00024407, -0.00029379], - [ 0.00018721, -0.00017669], - [ 0.00012562, -0.00009923], - [ 0.00005962, -0.00003655], - [ 0. , 0. ], - [ 0. , 0. ], - [ 0. , 0. ], - [ 0. , 0. ], - [ 0. , 0. ], + [0.0, 0.0], + [-0.00006594, -0.00004315], + [-0.00012749, -0.00009338], + [-0.00018649, -0.00018164], + [-0.00024325, -0.0002873], + [-0.00029591, -0.00043023], + [-0.00034592, -0.00058795], + [-0.00039183, -0.00077834], + [-0.00043521, -0.00098149], + [-0.00047474, -0.00121284], + [-0.00051191, -0.00145513], + [-0.0005455, -0.00172143], + [-0.0005769, -0.00199696], + [-0.00060501, -0.00229263], + [-0.00063112, -0.00259596], + [-0.00065425, -0.00291588], + [-0.00067557, -0.00324204], + [-0.00069422, -0.00358162], + [-0.00071126, -0.00392617], + [-0.00072592, -0.00428129], + [-0.00073916, -0.00464026], + [-0.00075031, -0.0050073], + [-0.00076022, -0.0053772], + [-0.00076831, -0.00575298], + [-0.00077534, -0.00613077], + [-0.00078084, -0.00651256], + [-0.00078546, -0.00689564], + [-0.00078884, -0.00728117], + [-0.00079152, -0.00766739], + [-0.00079327, -0.00805485], + [-0.00079451, -0.00844256], + [-0.00079513, -0.00883062], + [-0.00079547, -0.00921861], + [-0.0007956, -0.00960656], + [-0.00047553, -0.00960645], + [-0.00015576, -0.00960649], + [0.0001639, -0.00960666], + [0.00048371, -0.00960712], + [0.0008037, -0.00960768], + [0.00080342, -0.0092195], + [0.00080296, -0.00883138], + [0.0008021, -0.00844342], + [0.00080088, -0.00805553], + [0.00079888, -0.0076684], + [0.0007963, -0.00728171], + [0.00079264, -0.00689683], + [0.00078819, -0.00651292], + [0.00078235, -0.0061322], + [0.00077554, -0.00575311], + [0.00076705, -0.00537893], + [0.00075739, -0.00500716], + [0.00074577, -0.00464232], + [0.00073279, -0.00428083], + [0.00071756, -0.00392862], + [0.00070079, -0.00358081], + [0.00068145, -0.00324495], + [0.00066038, -0.00291468], + [0.00063644, -0.00259938], + [0.00061058, -0.00229101], + [0.00058154, -0.00200098], + [0.00055038, -0.00171937], + [0.00051574, -0.0014598], + [0.00047881, -0.0012103], + [0.00043813, -0.00098686], + [0.00039498, -0.00077525], + [0.00034784, -0.00059401], + [0.00029808, -0.00042646], + [0.00024407, -0.00029379], + [0.00018721, -0.00017669], + [0.00012562, -0.00009923], + [0.00005962, -0.00003655], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], ] test_suites.append((optimization_mtd_subg, expected_displacement_vector_subg)) @@ -249,7 +249,9 @@ def friction_bound(u_nu: float) -> float: setup = StaticSetup(mesh_descr) runner = StaticSolver(setup, solving_method) - result = runner.solve(initial_displacement=setup.initial_displacement, method=optimization_method) + result = runner.solve( + initial_displacement=setup.initial_displacement, method=optimization_method + ) displacement = result.body.mesh.nodes[:] - result.displaced_nodes[:] std_ids = standard_boundary_nodes(runner.body.mesh.nodes, runner.body.mesh.elements) From ac8587753d3125f553cd1599fa1a46b5863f401d Mon Sep 17 00:00:00 2001 From: Piotr Bartman-Szwarc Date: Fri, 26 Jul 2024 11:34:51 +0200 Subject: [PATCH 15/34] cleanup subgradient calc --- conmech/solvers/optimization/optimization.py | 10 ++++----- conmech/solvers/solver_methods.py | 23 +++++++------------- 2 files changed, 12 insertions(+), 21 deletions(-) diff --git a/conmech/solvers/optimization/optimization.py b/conmech/solvers/optimization/optimization.py index 76e4e857..8c3eed0d 100644 --- a/conmech/solvers/optimization/optimization.py +++ b/conmech/solvers/optimization/optimization.py @@ -34,7 +34,7 @@ from conmech.solvers.solver_methods import ( make_cost_functional, make_equation, - make_cost_functional_subgradient, + make_subgradient, ) @@ -62,11 +62,9 @@ def __init__( variable_dimension=statement.dimension_out, problem_dimension=statement.dimension_in, ) - if hasattr(contact_law, "subderivative_normal_direction"): - self.subgradient = make_cost_functional_subgradient( - djn=contact_law.subderivative_normal_direction, # TODO - djt=None, - dh_functional=None, + if hasattr(contact_law, "subderivative_normal_direction"): # TODO + self.subgradient = make_subgradient( + djn=contact_law.subderivative_normal_direction, ) if isinstance(statement, WaveStatement): if isinstance(contact_law, InteriorContactLaw): diff --git a/conmech/solvers/solver_methods.py b/conmech/solvers/solver_methods.py index f28f15f5..2fcafe77 100644 --- a/conmech/solvers/solver_methods.py +++ b/conmech/solvers/solver_methods.py @@ -161,9 +161,8 @@ def contact_cost(length, normal, normal_bound, tangential, tangential_bound): def contact_cost_functional( var, var_old, static_displacement, nodes, contact_boundary, contact_normals, dt ): - offset = len(var) // problem_dimension - cost = 0.0 + # pylint: disable=not-an-iterable for ei in numba.prange(len(contact_boundary)): edge = contact_boundary[ei] @@ -188,10 +187,6 @@ def contact_cost_functional( static_displacement_mean - static_displacement_normal * normal_vector ) - for node_id in edge: - if node_id >= offset: - continue - cost += contact_cost( nph.length(edge, nodes), normal_condition(vm_normal, static_displacement_normal, dt), @@ -226,18 +221,16 @@ def cost_functional( return cost_functional -def make_cost_functional_subgradient( - djn: Callable, djt: Optional[Callable] = None, dh_functional: Optional[Callable] = None +def make_subgradient( + djn: Callable, + problem_dimension=2, ): djn = njit(djn) - djt = njit(djt) - dh_functional = njit(dh_functional) - DIMENSION = 2 @numba.njit() - def contact_subgradient(u_vector, u_vector_old, nodes, contact_boundary, contact_normals): + def contact_subgradient(u_vector, nodes, contact_boundary): cost = np.zeros_like(u_vector) - offset = len(u_vector) // DIMENSION + offset = len(u_vector) // problem_dimension for edge in contact_boundary: n_id_0 = edge[0] @@ -254,7 +247,7 @@ def contact_subgradient(u_vector, u_vector_old, nodes, contact_boundary, contact cost[n_id_1 + offset] = cost[n_id_1] return cost - # pylint: disable=unused-argument # 'dt' + # pylint: disable=unused-argument # takes the same args as cost_functional @numba.njit() def subgradient( var, @@ -268,7 +261,7 @@ def subgradient( base_integrals, dt, ): - dj = contact_subgradient(var, var_old, nodes, contact_boundary, contact_normals) + dj = contact_subgradient(var, nodes, contact_boundary) result = np.dot(lhs, var) - rhs + dj result = result.ravel() return result From a0c2467e8102cbff2f0a46787f80f4ef45391577 Mon Sep 17 00:00:00 2001 From: Piotr Bartman-Szwarc Date: Wed, 31 Jul 2024 08:15:14 +0200 Subject: [PATCH 16/34] cleanup kosopt --- .../contact/damped_normal_compliance.py | 10 +++ conmech/plotting/membrane.py | 2 +- conmech/solvers/optimization/optimization.py | 25 ++++++- conmech/solvers/solver_methods.py | 8 ++- examples/Bagirrov_Bartman_Ochal_2024.py | 65 +++++++++---------- examples/Makela_et_al_1998.py | 2 +- .../regression/test_Makela_et_al_1998.py | 1 + 7 files changed, 69 insertions(+), 44 deletions(-) diff --git a/conmech/dynamics/contact/damped_normal_compliance.py b/conmech/dynamics/contact/damped_normal_compliance.py index de25c502..6b9d857e 100644 --- a/conmech/dynamics/contact/damped_normal_compliance.py +++ b/conmech/dynamics/contact/damped_normal_compliance.py @@ -52,4 +52,14 @@ def potential_normal_direction( return 0 return kappa * (displacement - obstacle_level) + beta * var_nu + @staticmethod + def subderivative_normal_direction( + var_nu: float, static_displacement_nu: float, dt: float + ) -> float: + displacement = static_displacement_nu + var_nu * dt + if displacement < obstacle_level: + return 0 + return kappa * dt + beta + + return DampedNormalCompliance diff --git a/conmech/plotting/membrane.py b/conmech/plotting/membrane.py index dd49f8ec..7111c3e3 100644 --- a/conmech/plotting/membrane.py +++ b/conmech/plotting/membrane.py @@ -45,7 +45,7 @@ def plot_limit_points( finish=True, ylim=(0, 1), ): - buff = np.zeros((2, 1000)) + buff = np.zeros((2, 1100)) # some extra space for multiple zeros buff_size = 0 for time, zeros in prod.data.items(): for zero in zeros: # significant performance improvement diff --git a/conmech/solvers/optimization/optimization.py b/conmech/solvers/optimization/optimization.py index 8c3eed0d..e37c4cd1 100644 --- a/conmech/solvers/optimization/optimization.py +++ b/conmech/solvers/optimization/optimization.py @@ -36,6 +36,7 @@ make_equation, make_subgradient, ) +from kosopt.qsmlm import make_minimizer class Optimization(Solver): @@ -66,6 +67,8 @@ def __init__( self.subgradient = make_subgradient( djn=contact_law.subderivative_normal_direction, ) + else: + self.subgradient = None if isinstance(statement, WaveStatement): if isinstance(contact_law, InteriorContactLaw): self.loss = make_equation( # TODO! @@ -81,6 +84,7 @@ def __init__( variable_dimension=statement.dimension_out, problem_dimension=statement.dimension_in, ) + self.minimizer = None def __str__(self): raise NotImplementedError() @@ -127,6 +131,16 @@ def _solve_impl( sols.append(solution) loss.append(self.loss(solution, *args)[0]) + if self.minimizer is None and method.lower() in ( + "quasi secant method", + "limited memory quasi secant method", + "quasi secant method limited memory", + "qsm", + "qsmlm", + "subgradient", + ): + self.minimizer = make_minimizer(self.loss, self.subgradient) + while norm >= fixed_point_abs_tol: if method.lower() in ( "quasi secant method", @@ -138,16 +152,20 @@ def _solve_impl( # pylint: disable=import-outside-toplevel,import-error) from kosopt import qsmlm - solution = qsmlm.minimize( - self.loss, solution, args=args, maxiter=maxiter, subgradient=self.subgradient + solution = self.minimizer( + solution, args, maxiter=maxiter ) + sols.append(solution.copy()) + loss.append(self.loss(solution, *args)[0]) elif method.lower() in ("subgradient",): # pylint: disable=import-outside-toplevel,import-error) from kosopt import subgradient solution = subgradient.minimize( - self.loss, solution, args=args, maxiter=maxiter, subgradient=self.subgradient + self.minimizer, self.loss, solution, args, maxiter=maxiter ) + sols.append(solution.copy()) + loss.append(self.loss(solution, *args)[0]) elif method.lower() in ( # TODO "discontinuous gradient", "discontinuous gradient method", @@ -158,6 +176,7 @@ def _solve_impl( solution = qsmlmi.minimize(self.loss, solution, args=args, maxiter=maxiter) sols.append(solution.copy()) + loss.append(self.loss(solution, *args)[0]) elif method.lower() == "constrained": contact_nodes_count = self.body.mesh.boundaries.contact_nodes_count diff --git a/conmech/solvers/solver_methods.py b/conmech/solvers/solver_methods.py index 2fcafe77..b9a074c8 100644 --- a/conmech/solvers/solver_methods.py +++ b/conmech/solvers/solver_methods.py @@ -261,9 +261,11 @@ def subgradient( base_integrals, dt, ): - dj = contact_subgradient(var, nodes, contact_boundary) - result = np.dot(lhs, var) - rhs + dj - result = result.ravel() + result = np.zeros_like(var) + ind = lhs.shape[0] + dj = contact_subgradient(var[:ind], nodes, contact_boundary) + result_ = np.dot(lhs, var[:ind]) - rhs + dj + result[:ind] = result_.ravel() return result return subgradient diff --git a/examples/Bagirrov_Bartman_Ochal_2024.py b/examples/Bagirrov_Bartman_Ochal_2024.py index 916dbafd..7d644c07 100644 --- a/examples/Bagirrov_Bartman_Ochal_2024.py +++ b/examples/Bagirrov_Bartman_Ochal_2024.py @@ -29,47 +29,40 @@ k30 = 2.5e12 * kN * surface -def normal_direction(u_nu: float) -> float: - if u_nu <= 0: - return 0.0 - if u_nu < 0.5 * mm: - return k0 * u_nu * 2 - if u_nu < 1 * mm: - return k10 * (u_nu * 2) + k11 - if u_nu < 2 * mm: - return k20 * (u_nu * 2) + k21 - return u_nu**3 * 4 * k30 - - class MMLV99(ContactLaw): @staticmethod - def potential_normal_direction(u_nu: float) -> float: - if u_nu <= 0: + def potential_normal_direction( + var_nu: float, static_displacement_nu: float, dt: float + ) -> float: + if var_nu <= 0: return 0.0 - if u_nu < 0.5 * mm: - return k0 * u_nu**2 - if u_nu < 1 * mm: - return k10 * u_nu**2 + k11 * u_nu - if u_nu < 2 * mm: - return k20 * u_nu**2 + k21 * u_nu + 4 - return u_nu**4 * k30 + if var_nu < 0.5 * mm: + return k0 * var_nu**2 + if var_nu < 1 * mm: + return k10 * var_nu**2 + k11 * var_nu + if var_nu < 2 * mm: + return k20 * var_nu**2 + k21 * var_nu + 4 + return var_nu**4 * k30 @staticmethod - def potential_tangential_direction(u_tau: np.ndarray) -> float: - return np.log(np.sum(u_tau * u_tau) ** 0.5 + 1) - - @staticmethod - def subderivative_normal_direction(u_nu: float, v_nu: float) -> float: - return 0 + def potential_tangential_direction( + var_tau: float, static_displacement_tau: float, dt: float + ) -> float: + return np.log(np.sum(var_tau ** 2) ** 0.5 + 1) @staticmethod - def regularized_subderivative_tangential_direction( - u_tau: np.ndarray, v_tau: np.ndarray, rho=1e-7 + def subderivative_normal_direction( + var_nu: float, static_displacement_nu: float, dt: float ) -> float: - """ - Coulomb regularization - """ - return 0 + if var_nu <= 0: + return 0.0 + if var_nu < 0.5 * mm: + return k0 * var_nu * 2 + if var_nu < 1 * mm: + return k10 * (var_nu * 2) + k11 + if var_nu < 2 * mm: + return k20 * (var_nu * 2) + k21 + return var_nu ** 3 * 4 * k30 @dataclass() @@ -173,7 +166,7 @@ def outer_forces(x, t=None): ) x = state.body.mesh.nodes[: state.body.mesh.contact_nodes_count - 1, 0] u = state.displacement[: state.body.mesh.contact_nodes_count - 1, 1] - y1 = [normal_direction(-u_) for u_ in u] + y1 = [MMLV99().subderivative_normal_direction(-u_, 0., 0.) for u_ in u] print(f) plt.plot(x, y1, label=f"{f:.2e}") plt.title(m) @@ -185,11 +178,11 @@ def outer_forces(x, t=None): X = np.linspace((2 - 2) * mm, (2 + 2) * mm, 1000) Y = np.empty(1000) for i in range(1000): - Y[i] = MMLV99.potential_normal_direction(X[i]) + Y[i] = MMLV99.potential_normal_direction(X[i], 0., 0.) plt.plot(X, Y) plt.show() for i in range(1000): - Y[i] = normal_direction(X[i]) + Y[i] = MMLV99().subderivative_normal_direction(X[i], 0., 0.) plt.plot(X, Y) plt.show() # results = { diff --git a/examples/Makela_et_al_1998.py b/examples/Makela_et_al_1998.py index 070f0998..caae931e 100644 --- a/examples/Makela_et_al_1998.py +++ b/examples/Makela_et_al_1998.py @@ -183,7 +183,7 @@ def outer_forces(x, t=None): ) x = state.body.mesh.nodes[: state.body.mesh.contact_nodes_count - 1, 0] u = state.displacement[: state.body.mesh.contact_nodes_count - 1, 1] - y1 = [MMLV99().normal_direction(-u_) for u_ in u] + y1 = [MMLV99().subderivative_normal_direction(-u_, 0.0, 0.0) for u_ in u] plt.plot(x, y1, label=f"{f:.2e}") plt.ylabel("Interlaminar binding force [kN/m$^2$]") plt.xlabel(r"Contact interface [mm]") diff --git a/tests/test_conmech/regression/test_Makela_et_al_1998.py b/tests/test_conmech/regression/test_Makela_et_al_1998.py index e6d29e1c..3a572508 100644 --- a/tests/test_conmech/regression/test_Makela_et_al_1998.py +++ b/tests/test_conmech/regression/test_Makela_et_al_1998.py @@ -201,6 +201,7 @@ def generate_test_suits(): ] test_suites.append((optimization_mtd_subg, expected_displacement_vector_subg)) + test_suites.append(("qsmlm", expected_displacement_vector_subg)) return test_suites From 502e52d26b3fd69f53f2a4915aa5503b8c83300e Mon Sep 17 00:00:00 2001 From: Piotr Bartman-Szwarc Date: Wed, 31 Jul 2024 08:25:31 +0200 Subject: [PATCH 17/34] black --- .../contact/damped_normal_compliance.py | 3 +-- conmech/solvers/optimization/optimization.py | 19 +++++++------------ examples/Bagirrov_Bartman_Ochal_2024.py | 10 +++++----- 3 files changed, 13 insertions(+), 19 deletions(-) diff --git a/conmech/dynamics/contact/damped_normal_compliance.py b/conmech/dynamics/contact/damped_normal_compliance.py index 6b9d857e..eebc378f 100644 --- a/conmech/dynamics/contact/damped_normal_compliance.py +++ b/conmech/dynamics/contact/damped_normal_compliance.py @@ -54,12 +54,11 @@ def potential_normal_direction( @staticmethod def subderivative_normal_direction( - var_nu: float, static_displacement_nu: float, dt: float + var_nu: float, static_displacement_nu: float, dt: float ) -> float: displacement = static_displacement_nu + var_nu * dt if displacement < obstacle_level: return 0 return kappa * dt + beta - return DampedNormalCompliance diff --git a/conmech/solvers/optimization/optimization.py b/conmech/solvers/optimization/optimization.py index e37c4cd1..a0af69a7 100644 --- a/conmech/solvers/optimization/optimization.py +++ b/conmech/solvers/optimization/optimization.py @@ -132,12 +132,12 @@ def _solve_impl( loss.append(self.loss(solution, *args)[0]) if self.minimizer is None and method.lower() in ( - "quasi secant method", - "limited memory quasi secant method", - "quasi secant method limited memory", - "qsm", - "qsmlm", - "subgradient", + "quasi secant method", + "limited memory quasi secant method", + "quasi secant method limited memory", + "qsm", + "qsmlm", + "subgradient", ): self.minimizer = make_minimizer(self.loss, self.subgradient) @@ -149,12 +149,7 @@ def _solve_impl( "qsm", "qsmlm", ): - # pylint: disable=import-outside-toplevel,import-error) - from kosopt import qsmlm - - solution = self.minimizer( - solution, args, maxiter=maxiter - ) + solution = self.minimizer(solution, args, maxiter=maxiter) sols.append(solution.copy()) loss.append(self.loss(solution, *args)[0]) elif method.lower() in ("subgradient",): diff --git a/examples/Bagirrov_Bartman_Ochal_2024.py b/examples/Bagirrov_Bartman_Ochal_2024.py index 7d644c07..3dd3534a 100644 --- a/examples/Bagirrov_Bartman_Ochal_2024.py +++ b/examples/Bagirrov_Bartman_Ochal_2024.py @@ -48,7 +48,7 @@ def potential_normal_direction( def potential_tangential_direction( var_tau: float, static_displacement_tau: float, dt: float ) -> float: - return np.log(np.sum(var_tau ** 2) ** 0.5 + 1) + return np.log(np.sum(var_tau**2) ** 0.5 + 1) @staticmethod def subderivative_normal_direction( @@ -62,7 +62,7 @@ def subderivative_normal_direction( return k10 * (var_nu * 2) + k11 if var_nu < 2 * mm: return k20 * (var_nu * 2) + k21 - return var_nu ** 3 * 4 * k30 + return var_nu**3 * 4 * k30 @dataclass() @@ -166,7 +166,7 @@ def outer_forces(x, t=None): ) x = state.body.mesh.nodes[: state.body.mesh.contact_nodes_count - 1, 0] u = state.displacement[: state.body.mesh.contact_nodes_count - 1, 1] - y1 = [MMLV99().subderivative_normal_direction(-u_, 0., 0.) for u_ in u] + y1 = [MMLV99().subderivative_normal_direction(-u_, 0.0, 0.0) for u_ in u] print(f) plt.plot(x, y1, label=f"{f:.2e}") plt.title(m) @@ -178,11 +178,11 @@ def outer_forces(x, t=None): X = np.linspace((2 - 2) * mm, (2 + 2) * mm, 1000) Y = np.empty(1000) for i in range(1000): - Y[i] = MMLV99.potential_normal_direction(X[i], 0., 0.) + Y[i] = MMLV99.potential_normal_direction(X[i], 0.0, 0.0) plt.plot(X, Y) plt.show() for i in range(1000): - Y[i] = MMLV99().subderivative_normal_direction(X[i], 0., 0.) + Y[i] = MMLV99().subderivative_normal_direction(X[i], 0.0, 0.0) plt.plot(X, Y) plt.show() # results = { From 7daf8f643e743574a5672fe3c72cabdccf2fb960 Mon Sep 17 00:00:00 2001 From: Piotr Bartman-Szwarc Date: Wed, 28 Aug 2024 10:20:28 +0200 Subject: [PATCH 18/34] WIP --- conmech/solvers/optimization/optimization.py | 1 + tests/test_conmech/regression/test_static.py | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/conmech/solvers/optimization/optimization.py b/conmech/solvers/optimization/optimization.py index a0af69a7..86c68868 100644 --- a/conmech/solvers/optimization/optimization.py +++ b/conmech/solvers/optimization/optimization.py @@ -212,6 +212,7 @@ def constr(x): norm = np.linalg.norm(np.subtract(solution, old_solution)) old_solution = solution.copy() min_index = loss.index(np.min(loss)) + print(np.min(loss)) # TODO solution = sols[min_index] return solution diff --git a/tests/test_conmech/regression/test_static.py b/tests/test_conmech/regression/test_static.py index 2efe6f60..bd36e4bc 100644 --- a/tests/test_conmech/regression/test_static.py +++ b/tests/test_conmech/regression/test_static.py @@ -15,7 +15,7 @@ from tests.test_conmech.regression.std_boundary import standard_boundary_nodes -@pytest.fixture(params=["direct", "global optimization", "schur"]) +@pytest.fixture(params=["direct", "global optimization", "schur"][-1:]) def solving_method(request): return request.param @@ -161,7 +161,7 @@ def outer_forces(x, t=None): @pytest.mark.parametrize("setup, expected_displacement_vector", generate_test_suits()) def test_static_solver(solving_method, setup, expected_displacement_vector): runner = StaticSolver(setup, solving_method) - result = runner.solve(initial_displacement=setup.initial_displacement) + result = runner.solve(initial_displacement=setup.initial_displacement, method="qsm") displacement = result.body.mesh.nodes[:] - result.displaced_nodes[:] std_ids = standard_boundary_nodes(runner.body.mesh.nodes, runner.body.mesh.elements) From 0786538d0fc720e6847ba2ec21fa2162368f9efb Mon Sep 17 00:00:00 2001 From: Piotr Bartman-Szwarc Date: Thu, 12 Sep 2024 15:52:41 +0200 Subject: [PATCH 19/34] fix tests make schur respect fixed_point_abs_tol param full subgradient computation --- .../contact/relu_slope_contact_law.py | 2 +- conmech/solvers/optimization/optimization.py | 7 +- .../solvers/optimization/schur_complement.py | 2 - conmech/solvers/solver_methods.py | 94 +++++++++++++------ .../test_conmech/regression/test_static_3d.py | 52 +++++----- 5 files changed, 98 insertions(+), 59 deletions(-) diff --git a/conmech/dynamics/contact/relu_slope_contact_law.py b/conmech/dynamics/contact/relu_slope_contact_law.py index 7e2e8513..f77a897c 100644 --- a/conmech/dynamics/contact/relu_slope_contact_law.py +++ b/conmech/dynamics/contact/relu_slope_contact_law.py @@ -27,6 +27,6 @@ def subderivative_normal_direction( ) -> float: if var_nu <= 0: return 0.0 - return slope * var_nu * static_displacement_nu + return slope * var_nu return ReLuContactLaw diff --git a/conmech/solvers/optimization/optimization.py b/conmech/solvers/optimization/optimization.py index 86c68868..672f13c9 100644 --- a/conmech/solvers/optimization/optimization.py +++ b/conmech/solvers/optimization/optimization.py @@ -65,7 +65,7 @@ def __init__( ) if hasattr(contact_law, "subderivative_normal_direction"): # TODO self.subgradient = make_subgradient( - djn=contact_law.subderivative_normal_direction, + normal_condition=contact_law.subderivative_normal_direction, ) else: self.subgradient = None @@ -157,7 +157,8 @@ def _solve_impl( from kosopt import subgradient solution = subgradient.minimize( - self.minimizer, self.loss, solution, args, maxiter=maxiter + self.minimizer, self.loss, solution, args, + maxiter=maxiter, subgradient=self.subgradient ) sols.append(solution.copy()) loss.append(self.loss(solution, *args)[0]) @@ -212,7 +213,7 @@ def constr(x): norm = np.linalg.norm(np.subtract(solution, old_solution)) old_solution = solution.copy() min_index = loss.index(np.min(loss)) - print(np.min(loss)) # TODO + print(method, np.min(loss)) # TODO solution = sols[min_index] return solution diff --git a/conmech/solvers/optimization/schur_complement.py b/conmech/solvers/optimization/schur_complement.py index 58658815..cfbcff60 100644 --- a/conmech/solvers/optimization/schur_complement.py +++ b/conmech/solvers/optimization/schur_complement.py @@ -77,8 +77,6 @@ def rhs(self) -> np.ndarray: def _solve_impl( self, initial_guess: np.ndarray, - *, - fixed_point_abs_tol: float = math.inf, **kwargs, ) -> np.ndarray: # initial_guess: [[xc, xf, xd] [yc, yf, yd]] diff --git a/conmech/solvers/solver_methods.py b/conmech/solvers/solver_methods.py index b9a074c8..2d12aff0 100644 --- a/conmech/solvers/solver_methods.py +++ b/conmech/solvers/solver_methods.py @@ -32,10 +32,10 @@ def interpolate_node_between(edge, vector, full_vector, dimension): for node in edge: for i in range(dimension): if node < offset: # exclude dirichlet nodes (and inner nodes in schur) - result[i] += 0.5 * vector[i * offset + node] + result[i] += vector[i * offset + node] / len(edge) else: # old values - result[i] += 0.5 * full_vector[i * offset_full + node] + result[i] += full_vector[i * offset_full + node] / len(edge) return result @@ -94,8 +94,8 @@ def contact_part(u_vector, nodes, contact_boundary, contact_normals): um_normal = (um * normal_vector).sum() edge_len = nph.length(edge, nodes) - j_x = edge_len * 0.5 * (jn(um_normal, normal_vector[0], 0.0)) - j_y = edge_len * 0.5 * (jn(um_normal, normal_vector[1], 0.0)) + j_x = edge_len * 0.5 * (jn(um_normal, 0.0, 0.0) * normal_vector[0]) + j_y = edge_len * 0.5 * (jn(um_normal, 0.0, 0.0) * normal_vector[1]) if n_id_0 < offset: contact_vector[n_id_0] += j_x @@ -222,34 +222,72 @@ def cost_functional( def make_subgradient( - djn: Callable, + normal_condition: Callable, + normal_condition_bound: Optional[Callable] = None, + tangential_condition: Optional[Callable] = None, + tangential_condition_bound: Optional[Callable] = None, problem_dimension=2, + variable_dimension=2, ): - djn = njit(djn) + normal_condition = njit(normal_condition) + normal_condition_bound = njit(normal_condition_bound, value=1) + tangential_condition = njit(tangential_condition) + tangential_condition_bound = njit(tangential_condition_bound, value=1) @numba.njit() - def contact_subgradient(u_vector, nodes, contact_boundary): - cost = np.zeros_like(u_vector) - offset = len(u_vector) // problem_dimension - - for edge in contact_boundary: - n_id_0 = edge[0] - n_id_1 = edge[1] - n_0 = nodes[n_id_0] - n_1 = nodes[n_id_1] - if n_id_0 < offset: - um_normal_0 = -n_0[0] # TODO - cost[n_id_0] = djn(um_normal_0, 0.0, 0.0) - cost[n_id_0 + offset] = cost[n_id_0] - if n_id_1 < offset: - um_normal_1 = -n_1[0] # TODO - cost[n_id_1] = djn(um_normal_1, 0.0, 0.0) - cost[n_id_1 + offset] = cost[n_id_1] + def contact_cost(length, normal, normal_bound, tangential, tangential_bound): + return length * (normal_bound * normal + tangential_bound * tangential) + + @numba.njit() + def contact_subgradient( + var, var_old, static_displacement, nodes, contact_boundary, contact_normals, dt + ): + cost = np.zeros_like(var) + offset = len(var) // variable_dimension + + # pylint: disable=not-an-iterable + for ei in numba.prange(len(contact_boundary)): + edge = contact_boundary[ei] + normal_vector = contact_normals[ei] + # ASSUMING `u_vector` and `nodes` have the same order! + vm = interpolate_node_between(edge, var, var_old, dimension=variable_dimension) + if variable_dimension == 1: + raise NotImplementedError() + vm_normal = vm[0] + vm_tangential = np.empty(0) + else: + vm_normal = (vm * normal_vector).sum() + vm_tangential = vm - vm_normal * normal_vector + + static_displacement_mean = interpolate_node_between( + edge, + static_displacement, + static_displacement, + dimension=problem_dimension, + ) + static_displacement_normal = (static_displacement_mean * normal_vector).sum() + static_displacement_tangential = ( + static_displacement_mean - static_displacement_normal * normal_vector + ) + + subgrad = contact_cost( + nph.length(edge, nodes), + normal_condition(vm_normal, static_displacement_normal, dt), + normal_condition_bound(vm_normal, static_displacement_normal, dt), + tangential_condition(vm_tangential, static_displacement_tangential, dt), + tangential_condition_bound(vm_normal, static_displacement_normal, dt), + ) + + for node in edge: + for i in range(variable_dimension): + if node < offset: + cost[i * offset + node] += normal_vector[i] / len(edge) * subgrad + return cost - # pylint: disable=unused-argument # takes the same args as cost_functional + # pylint: disable=too-many-arguments,unused-argument # 'base_integrals' @numba.njit() - def subgradient( + def cost_functional( var, var_old, nodes, @@ -262,10 +300,12 @@ def subgradient( dt, ): result = np.zeros_like(var) + dj = contact_subgradient( + var, var_old, u_vector, nodes, contact_boundary, contact_normals, dt + ) ind = lhs.shape[0] - dj = contact_subgradient(var[:ind], nodes, contact_boundary) result_ = np.dot(lhs, var[:ind]) - rhs + dj result[:ind] = result_.ravel() return result - return subgradient + return cost_functional diff --git a/tests/test_conmech/regression/test_static_3d.py b/tests/test_conmech/regression/test_static_3d.py index 69b82909..afc48e60 100644 --- a/tests/test_conmech/regression/test_static_3d.py +++ b/tests/test_conmech/regression/test_static_3d.py @@ -50,32 +50,32 @@ def outer_forces(x, t=None): setup_m02_m02 = StaticSetup(mesh_descr) expected_displacement_vector_m02_m02 = [ - [0.0, 0.0, 0.0], - [0.0, 0.0, 0.0], - [0.0, 0.0, 0.0], - [0.00886805, 0.00759973, 0.01309122], - [0.01251517, 0.01304646, 0.02180133], - [0.0047388, 0.01377955, 0.02134457], - [-0.00179084, 0.01481824, 0.02136016], - [0.00211174, 0.01387662, 0.02174084], - [0.00165348, 0.00748375, 0.01153221], - [-0.00196903, 0.00853421, 0.01211016], - [0.0, 0.0, 0.0], - [0.0, 0.0, 0.0], - [-0.00548426, 0.00962111, 0.01202977], - [-0.0060029, 0.01529378, 0.02094435], - [0.00087541, 0.01483005, 0.02075665], - [-0.00010425, 0.00850936, 0.01158434], - [0.0, 0.0, 0.0], - [0.0, 0.0, 0.0], - [0.00577241, 0.00703007, 0.01156364], - [0.00711386, 0.0132022, 0.02045272], - [0.0, 0.0, 0.0], - [0.0, 0.0, 0.0], - [0.01268809, 0.0081164, 0.0137591], - [0.01424769, 0.01004692, 0.02280081], - [0.00928219, 0.01244232, 0.02200091], - [0.00655036, 0.00806369, 0.01241398], + [0., 0., 0.], + [0., 0., 0.], + [0., 0., 0.], + [0.00887781, 0.00938318, 0.01303338], + [0.01254745, 0.0164677, 0.02173087], + [0.00484257, 0.01694291, 0.02137632], + [-0.00172692, 0.01755931, 0.02152949], + [0.00322087, 0.01678243, 0.0214544], + [0.00253174, 0.00921583, 0.01132137], + [-0.00196061, 0.00997206, 0.01218131], + [0., 0., 0.], + [0., 0., 0.], + [-0.0062681, 0.01092748, 0.0123151], + [-0.00694282, 0.01796747, 0.021246], + [-0.00016377, 0.01769789, 0.02100172], + [-0.00086678, 0.0100699, 0.0117306], + [0., 0., 0.], + [0., 0., 0.], + [0.00494897, 0.00864805, 0.0115685], + [0.00615472, 0.0164003, 0.02059693], + [0., 0., 0.], + [0., 0., 0.], + [0.01361542, 0.01025, 0.01362368], + [0.01573552, 0.01474276, 0.02232607], + [0.01065649, 0.01625437, 0.02177309], + [0.00752115, 0.01009524, 0.0122507], ] test_suites.append((setup_m02_m02, expected_displacement_vector_m02_m02)) From 356613a729e15bd59b61ebde79951cd5eb2a7cf1 Mon Sep 17 00:00:00 2001 From: Piotr Bartman-Szwarc Date: Thu, 12 Sep 2024 16:52:53 +0200 Subject: [PATCH 20/34] update tests --- tests/test_conmech/regression/test_static.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/tests/test_conmech/regression/test_static.py b/tests/test_conmech/regression/test_static.py index bd36e4bc..52f056c5 100644 --- a/tests/test_conmech/regression/test_static.py +++ b/tests/test_conmech/regression/test_static.py @@ -15,11 +15,16 @@ from tests.test_conmech.regression.std_boundary import standard_boundary_nodes -@pytest.fixture(params=["direct", "global optimization", "schur"][-1:]) +@pytest.fixture(params=["direct", "global optimization", "schur"]) def solving_method(request): return request.param +@pytest.fixture(params=["BFGS", "qsm", "subgradient"]) +def opt_method(request): + return request.param + + def generate_test_suits(): test_suites = [] @@ -159,9 +164,9 @@ def outer_forces(x, t=None): @pytest.mark.parametrize("setup, expected_displacement_vector", generate_test_suits()) -def test_static_solver(solving_method, setup, expected_displacement_vector): +def test_static_solver(solving_method, opt_method, setup, expected_displacement_vector): runner = StaticSolver(setup, solving_method) - result = runner.solve(initial_displacement=setup.initial_displacement, method="qsm") + result = runner.solve(initial_displacement=setup.initial_displacement, method=opt_method) displacement = result.body.mesh.nodes[:] - result.displaced_nodes[:] std_ids = standard_boundary_nodes(runner.body.mesh.nodes, runner.body.mesh.elements) From bf891135e6177954f95cf9c9d8a4ae77517f3274 Mon Sep 17 00:00:00 2001 From: Piotr Bartman-Szwarc Date: Thu, 12 Sep 2024 16:55:18 +0200 Subject: [PATCH 21/34] update Makela et al contact law --- examples/Makela_et_al_1998.py | 33 +++++++++++++++++++++------------ 1 file changed, 21 insertions(+), 12 deletions(-) diff --git a/examples/Makela_et_al_1998.py b/examples/Makela_et_al_1998.py index caae931e..ee3fc52c 100644 --- a/examples/Makela_et_al_1998.py +++ b/examples/Makela_et_al_1998.py @@ -49,31 +49,30 @@ def potential_normal_direction( var_nu: float, static_displacement_nu: float, dt: float ) -> float: u_nu = -var_nu - coef = 1.0 - if u_nu <= 0: + if u_nu <= 0.0: return 0.0 if u_nu < 0.5 * mm: - return k0 * u_nu**2 * coef - if u_nu < 1 * mm: - return (k10 * u_nu**2 + k11 * u_nu) * coef - if u_nu < 2 * mm: - return (k20 * u_nu**2 + k21 * u_nu + 4) * coef - return 16 * coef + return k0 * u_nu**2 + if u_nu < 1.0 * mm: + return k10 * u_nu**2 + k11 * u_nu + if u_nu < 2.0 * mm: + return k20 * u_nu**2 + k21 * u_nu + 4 + return 16.0 @staticmethod def subderivative_normal_direction( var_nu: float, static_displacement_nu: float, dt: float ) -> float: u_nu = -var_nu - if u_nu <= 0: + if u_nu <= 0.0: return 0.0 if u_nu < 0.5 * mm: return k0 * u_nu * 2 - if u_nu < 1 * mm: + if u_nu < 1.0 * mm: return k10 * (u_nu * 2) + k11 - if u_nu < 2 * mm: + if u_nu < 2.0 * mm: return k20 * (u_nu * 2) + k21 - return 0 + return 0.0 @staticmethod def potential_tangential_direction( @@ -81,6 +80,16 @@ def potential_tangential_direction( ) -> float: return np.log(np.sum(var_tau**2) ** 0.5 + 1) + @staticmethod + def subderivative_tangential_direction( + var_tau: float, static_displacement_tau: float, dt: float + ) -> float: + quadsum = np.sum(var_tau ** 2) + norm = quadsum ** 0.5 + denom = norm + quadsum + coef = 1 / denom if denom != 0.0 else 0.0 + return var_tau * coef + @dataclass() class StaticSetup(StaticDisplacementProblem): From 6c3d5c471a77da4d48de948c8ba06fd729521292 Mon Sep 17 00:00:00 2001 From: Piotr Bartman-Szwarc Date: Thu, 12 Sep 2024 17:20:00 +0200 Subject: [PATCH 22/34] small cleanup --- conmech/solvers/optimization/optimization.py | 47 ++++++++++--------- .../solvers/optimization/schur_complement.py | 1 - conmech/solvers/solver_methods.py | 12 ++--- examples/Bagirrov_Bartman_Ochal_2024.py | 2 +- examples/Makela_et_al_1998.py | 8 ++-- .../regression/test_Makela_et_al_1998.py | 2 +- tests/test_conmech/regression/test_static.py | 2 +- .../test_conmech/regression/test_static_3d.py | 18 +++---- 8 files changed, 48 insertions(+), 44 deletions(-) diff --git a/conmech/solvers/optimization/optimization.py b/conmech/solvers/optimization/optimization.py index 672f13c9..68d947e5 100644 --- a/conmech/solvers/optimization/optimization.py +++ b/conmech/solvers/optimization/optimization.py @@ -39,6 +39,22 @@ from kosopt.qsmlm import make_minimizer +QSMLM_NAMES = { + "quasi secant method", + "limited memory quasi secant method", + "quasi secant method limited memory", + "qsm", + "qsmlm", +} +GLOBAL_QSMLM_NAMES = { + "global quasi secant method", + "global limited memory quasi secant method", + "global quasi secant method limited memory", + "globqsm", + "globqsmlm", +} + + class Optimization(Solver): def __init__( self, @@ -131,34 +147,25 @@ def _solve_impl( sols.append(solution) loss.append(self.loss(solution, *args)[0]) - if self.minimizer is None and method.lower() in ( - "quasi secant method", - "limited memory quasi secant method", - "quasi secant method limited memory", - "qsm", - "qsmlm", - "subgradient", - ): + if self.minimizer is None and method.lower() in QSMLM_NAMES.union(GLOBAL_QSMLM_NAMES): self.minimizer = make_minimizer(self.loss, self.subgradient) while norm >= fixed_point_abs_tol: - if method.lower() in ( - "quasi secant method", - "limited memory quasi secant method", - "quasi secant method limited memory", - "qsm", - "qsmlm", - ): + if method.lower() in QSMLM_NAMES: solution = self.minimizer(solution, args, maxiter=maxiter) sols.append(solution.copy()) loss.append(self.loss(solution, *args)[0]) - elif method.lower() in ("subgradient",): + elif method.lower() in GLOBAL_QSMLM_NAMES: # pylint: disable=import-outside-toplevel,import-error) from kosopt import subgradient solution = subgradient.minimize( - self.minimizer, self.loss, solution, args, - maxiter=maxiter, subgradient=self.subgradient + self.minimizer, + self.loss, + solution, + args, + maxiter=maxiter, + subgradient=self.subgradient, ) sols.append(solution.copy()) loss.append(self.loss(solution, *args)[0]) @@ -195,7 +202,6 @@ def constr(x): solution = result.x sols.append(solution) loss.append(self.loss(solution, *args)[0]) - break else: result = scipy.optimize.minimize( self.loss, @@ -208,12 +214,11 @@ def constr(x): solution = result.x sols.append(solution.copy()) loss.append(self.loss(solution, *args)[0]) - break norm = np.linalg.norm(np.subtract(solution, old_solution)) old_solution = solution.copy() min_index = loss.index(np.min(loss)) - print(method, np.min(loss)) # TODO + print(method, np.min(loss)) # TODO solution = sols[min_index] return solution diff --git a/conmech/solvers/optimization/schur_complement.py b/conmech/solvers/optimization/schur_complement.py index cfbcff60..8475e7fa 100644 --- a/conmech/solvers/optimization/schur_complement.py +++ b/conmech/solvers/optimization/schur_complement.py @@ -2,7 +2,6 @@ Created at 22.02.2021 """ -import math from typing import Tuple import numpy as np diff --git a/conmech/solvers/solver_methods.py b/conmech/solvers/solver_methods.py index 2d12aff0..dfb2e256 100644 --- a/conmech/solvers/solver_methods.py +++ b/conmech/solvers/solver_methods.py @@ -252,12 +252,12 @@ def contact_subgradient( # ASSUMING `u_vector` and `nodes` have the same order! vm = interpolate_node_between(edge, var, var_old, dimension=variable_dimension) if variable_dimension == 1: - raise NotImplementedError() - vm_normal = vm[0] - vm_tangential = np.empty(0) - else: - vm_normal = (vm * normal_vector).sum() - vm_tangential = vm - vm_normal * normal_vector + raise NotImplementedError() # TODO + # vm_normal = vm[0] + # vm_tangential = np.empty(0) + # else: + vm_normal = (vm * normal_vector).sum() + vm_tangential = vm - vm_normal * normal_vector static_displacement_mean = interpolate_node_between( edge, diff --git a/examples/Bagirrov_Bartman_Ochal_2024.py b/examples/Bagirrov_Bartman_Ochal_2024.py index 3dd3534a..cb8378cb 100644 --- a/examples/Bagirrov_Bartman_Ochal_2024.py +++ b/examples/Bagirrov_Bartman_Ochal_2024.py @@ -228,7 +228,7 @@ def outer_forces(x, t=None): # plt.legend() # # plt.loglog() # plt.show() - methods = ("BFGS", "CG", "qsm", "Powell", "subgradient") + methods = ("BFGS", "CG", "qsm", "Powell", "globqsm") forces = ( 23e3 * kN, 25e3 * kN, diff --git a/examples/Makela_et_al_1998.py b/examples/Makela_et_al_1998.py index ee3fc52c..fe65f58e 100644 --- a/examples/Makela_et_al_1998.py +++ b/examples/Makela_et_al_1998.py @@ -82,10 +82,10 @@ def potential_tangential_direction( @staticmethod def subderivative_tangential_direction( - var_tau: float, static_displacement_tau: float, dt: float + var_tau: float, static_displacement_tau: float, dt: float ) -> float: - quadsum = np.sum(var_tau ** 2) - norm = quadsum ** 0.5 + quadsum = np.sum(var_tau**2) + norm = quadsum**0.5 denom = norm + quadsum coef = 1 / denom if denom != 0.0 else 0.0 return var_tau * coef @@ -273,5 +273,5 @@ def outer_forces(x, t=None): # plt.xlabel(r"Load [kN/m$^2$]") # plt.grid() # plt.show() - methods = ("BFGS", "CG", "Powell", "subgradient")[-2:] + methods = ("BFGS", "CG", "Powell", "globqsm")[-2:] main(Config(save=False, show=True, force=True).init(), methods, forces) diff --git a/tests/test_conmech/regression/test_Makela_et_al_1998.py b/tests/test_conmech/regression/test_Makela_et_al_1998.py index 3a572508..7b699e3b 100644 --- a/tests/test_conmech/regression/test_Makela_et_al_1998.py +++ b/tests/test_conmech/regression/test_Makela_et_al_1998.py @@ -120,7 +120,7 @@ def generate_test_suits(): ] test_suites.append((optimization_mtd_pow, expected_displacement_vector_pow)) - optimization_mtd_subg = "subgradient" + optimization_mtd_subg = "globqsm" expected_displacement_vector_subg = [ [0.0, 0.0], [-0.00006594, -0.00004315], diff --git a/tests/test_conmech/regression/test_static.py b/tests/test_conmech/regression/test_static.py index 52f056c5..f38bed6b 100644 --- a/tests/test_conmech/regression/test_static.py +++ b/tests/test_conmech/regression/test_static.py @@ -20,7 +20,7 @@ def solving_method(request): return request.param -@pytest.fixture(params=["BFGS", "qsm", "subgradient"]) +@pytest.fixture(params=["BFGS", "qsm", "globqsm"]) def opt_method(request): return request.param diff --git a/tests/test_conmech/regression/test_static_3d.py b/tests/test_conmech/regression/test_static_3d.py index afc48e60..61ca2d06 100644 --- a/tests/test_conmech/regression/test_static_3d.py +++ b/tests/test_conmech/regression/test_static_3d.py @@ -50,9 +50,9 @@ def outer_forces(x, t=None): setup_m02_m02 = StaticSetup(mesh_descr) expected_displacement_vector_m02_m02 = [ - [0., 0., 0.], - [0., 0., 0.], - [0., 0., 0.], + [0.0, 0.0, 0.0], + [0.0, 0.0, 0.0], + [0.0, 0.0, 0.0], [0.00887781, 0.00938318, 0.01303338], [0.01254745, 0.0164677, 0.02173087], [0.00484257, 0.01694291, 0.02137632], @@ -60,18 +60,18 @@ def outer_forces(x, t=None): [0.00322087, 0.01678243, 0.0214544], [0.00253174, 0.00921583, 0.01132137], [-0.00196061, 0.00997206, 0.01218131], - [0., 0., 0.], - [0., 0., 0.], + [0.0, 0.0, 0.0], + [0.0, 0.0, 0.0], [-0.0062681, 0.01092748, 0.0123151], [-0.00694282, 0.01796747, 0.021246], [-0.00016377, 0.01769789, 0.02100172], [-0.00086678, 0.0100699, 0.0117306], - [0., 0., 0.], - [0., 0., 0.], + [0.0, 0.0, 0.0], + [0.0, 0.0, 0.0], [0.00494897, 0.00864805, 0.0115685], [0.00615472, 0.0164003, 0.02059693], - [0., 0., 0.], - [0., 0., 0.], + [0.0, 0.0, 0.0], + [0.0, 0.0, 0.0], [0.01361542, 0.01025, 0.01362368], [0.01573552, 0.01474276, 0.02232607], [0.01065649, 0.01625437, 0.02177309], From c3de2f2a36a6c743367e70dd0d59538173a47a9f Mon Sep 17 00:00:00 2001 From: Piotr Bartman-Szwarc Date: Thu, 12 Sep 2024 17:27:17 +0200 Subject: [PATCH 23/34] fic CI --- conmech/solvers/optimization/optimization.py | 4 +++- tests/test_conmech/regression/test_static.py | 9 ++++++++- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/conmech/solvers/optimization/optimization.py b/conmech/solvers/optimization/optimization.py index 68d947e5..aa39e122 100644 --- a/conmech/solvers/optimization/optimization.py +++ b/conmech/solvers/optimization/optimization.py @@ -36,7 +36,6 @@ make_equation, make_subgradient, ) -from kosopt.qsmlm import make_minimizer QSMLM_NAMES = { @@ -148,6 +147,9 @@ def _solve_impl( loss.append(self.loss(solution, *args)[0]) if self.minimizer is None and method.lower() in QSMLM_NAMES.union(GLOBAL_QSMLM_NAMES): + # pylint: disable=import-outside-toplevel,import-error) + from kosopt.qsmlm import make_minimizer + self.minimizer = make_minimizer(self.loss, self.subgradient) while norm >= fixed_point_abs_tol: diff --git a/tests/test_conmech/regression/test_static.py b/tests/test_conmech/regression/test_static.py index f38bed6b..96fe6cf4 100644 --- a/tests/test_conmech/regression/test_static.py +++ b/tests/test_conmech/regression/test_static.py @@ -14,13 +14,20 @@ from conmech.dynamics.contact.relu_slope_contact_law import make_slope_contact_law from tests.test_conmech.regression.std_boundary import standard_boundary_nodes +try: + import kosopt + + available_opt_mtds = ["BFGS", "qsm", "globqsm"] +except ImportError: + available_opt_mtds = ["BFGS"] + @pytest.fixture(params=["direct", "global optimization", "schur"]) def solving_method(request): return request.param -@pytest.fixture(params=["BFGS", "qsm", "globqsm"]) +@pytest.fixture(params=available_opt_mtds) def opt_method(request): return request.param From a2e5f665a8dd30cf9507bd4e5016c8ee2613cdc0 Mon Sep 17 00:00:00 2001 From: Piotr Bartman-Szwarc Date: Thu, 12 Sep 2024 17:53:36 +0200 Subject: [PATCH 24/34] fix CI, drop support for python 3.9, add support for python 3.12 --- .github/workflows/conmech_tests.yml | 2 +- tests/test_conmech/regression/std_boundary.py | 2 +- .../regression/test_Makela_et_al_1998.py | 19 ++++++++++++++----- 3 files changed, 16 insertions(+), 7 deletions(-) diff --git a/.github/workflows/conmech_tests.yml b/.github/workflows/conmech_tests.yml index 8a9bb02a..b4bcd797 100644 --- a/.github/workflows/conmech_tests.yml +++ b/.github/workflows/conmech_tests.yml @@ -15,7 +15,7 @@ jobs: strategy: matrix: platform: [ ubuntu-latest, macos-13, windows-latest ] - python-version: [ "3.9", "3.10", "3.11"] + python-version: [ "3.10", "3.11", "3.12"] runs-on: ${{ matrix.platform }} steps: - uses: actions/checkout@v2 diff --git a/tests/test_conmech/regression/std_boundary.py b/tests/test_conmech/regression/std_boundary.py index 1cf4d44f..fa23f67d 100644 --- a/tests/test_conmech/regression/std_boundary.py +++ b/tests/test_conmech/regression/std_boundary.py @@ -7,7 +7,7 @@ def standard_boundary_nodes(nodes, elements): - """ + r""" Return nodes indices counter-clockwise for standard body (rectangle) with first node in (0, 0). For body: diff --git a/tests/test_conmech/regression/test_Makela_et_al_1998.py b/tests/test_conmech/regression/test_Makela_et_al_1998.py index 7b699e3b..8009a6c0 100644 --- a/tests/test_conmech/regression/test_Makela_et_al_1998.py +++ b/tests/test_conmech/regression/test_Makela_et_al_1998.py @@ -24,11 +24,17 @@ from conmech.mesh.boundaries_description import BoundariesDescription from conmech.scenarios.problems import StaticDisplacementProblem from conmech.simulations.problem_solver import StaticSolver -from conmech.properties.mesh_description import CrossMeshDescription, RectangleMeshDescription -from conmech.solvers.optimization.optimization import Optimization +from conmech.properties.mesh_description import RectangleMeshDescription from examples.Makela_et_al_1998 import MMLV99 from tests.test_conmech.regression.std_boundary import standard_boundary_nodes +try: + import kosopt + + available_opt_mtds = ["Powell", "qsm", "globqsm"] +except ImportError: + available_opt_mtds = ["Powell"] + @pytest.fixture(params=["schur"]) def solving_method(request): @@ -118,7 +124,8 @@ def generate_test_suits(): [0.0, 0.0], [0.0, 0.0], ] - test_suites.append((optimization_mtd_pow, expected_displacement_vector_pow)) + if optimization_mtd_pow in available_opt_mtds: + test_suites.append((optimization_mtd_pow, expected_displacement_vector_pow)) optimization_mtd_subg = "globqsm" expected_displacement_vector_subg = [ @@ -200,8 +207,10 @@ def generate_test_suits(): [0.0, 0.0], ] - test_suites.append((optimization_mtd_subg, expected_displacement_vector_subg)) - test_suites.append(("qsmlm", expected_displacement_vector_subg)) + if optimization_mtd_subg in available_opt_mtds: + test_suites.append((optimization_mtd_subg, expected_displacement_vector_subg)) + if "qsmlm" in available_opt_mtds: + test_suites.append(("qsmlm", expected_displacement_vector_subg)) return test_suites From 89e005597f3c54d4cb345014d93f3cd7d3bdd69f Mon Sep 17 00:00:00 2001 From: Piotr Bartman-Szwarc Date: Thu, 12 Sep 2024 18:03:28 +0200 Subject: [PATCH 25/34] upgrade matplotlib to fix CI at windows --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index cc23f6f8..ef246f0e 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,7 +1,7 @@ numpy==2.0.1 numba==0.60.0 jupyter==1.0.0 -matplotlib==3.9.1 +matplotlib==3.9.2 scipy==1.13.1 pytest==8.3.2 pygmsh==7.1.17 From 2918075d0260c2827d9ddca63c4ba743e81431c2 Mon Sep 17 00:00:00 2001 From: Piotr Bartman-Szwarc Date: Thu, 12 Sep 2024 19:39:30 +0200 Subject: [PATCH 26/34] cleanup --- examples/Sofonea_Ochal_Bartman_2023.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/examples/Sofonea_Ochal_Bartman_2023.py b/examples/Sofonea_Ochal_Bartman_2023.py index 4d83a53a..306c139d 100644 --- a/examples/Sofonea_Ochal_Bartman_2023.py +++ b/examples/Sofonea_Ochal_Bartman_2023.py @@ -167,12 +167,6 @@ def zero_relaxation(t=None): setup.relaxation = examples[name]["relaxation"] runner = QuasistaticRelaxation(setup, solving_method="schur") - bid = runner.body.mesh.contact_boundary[:, 0] - eid = runner.body.mesh.contact_boundary[:, 1] - bx = runner.body.mesh.nodes[bid][:, 0] - ex = runner.body.mesh.nodes[eid][:, 0] - length = np.max(ex - bx) - print(length) states = runner.solve( n_steps=examples[name]["n_steps"], From 460b406cf6eb2a4869ca1c68b8d5e31db20ee5da Mon Sep 17 00:00:00 2001 From: Piotr Bartman-Szwarc Date: Thu, 12 Sep 2024 20:12:53 +0200 Subject: [PATCH 27/34] minor fix --- examples/BOST_2024.py | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/BOST_2024.py b/examples/BOST_2024.py index 23f483b5..222cd873 100644 --- a/examples/BOST_2024.py +++ b/examples/BOST_2024.py @@ -214,6 +214,7 @@ def main( def plot_errors(config, X, Y, highlighted_id, save: Optional[str] = None): plt.plot(X[:], Y[:], marker="o", color="gray") if highlighted_id is not None: + highlighted_id = list(highlighted_id) plt.plot(X[highlighted_id], Y[highlighted_id], "ro", color="black") plt.loglog() plt.grid(True, which="major", linestyle="--", linewidth=0.5) From ad1117973a9f4e6cea181d2dcab429dedb436ab0 Mon Sep 17 00:00:00 2001 From: Piotr Bartman-Szwarc Date: Tue, 24 Sep 2024 17:58:12 +0200 Subject: [PATCH 28/34] cleanup --- conmech/plotting/drawer.py | 23 +++- conmech/solvers/optimization/optimization.py | 6 +- examples/Bagirrov_Bartman_Ochal_2024.py | 69 ++++------ examples/Makela_et_al_1998.py | 125 +++++++++++-------- 4 files changed, 114 insertions(+), 109 deletions(-) diff --git a/conmech/plotting/drawer.py b/conmech/plotting/drawer.py index 43af9857..d2d3f132 100644 --- a/conmech/plotting/drawer.py +++ b/conmech/plotting/drawer.py @@ -55,6 +55,22 @@ def draw( save_format="png", title=None, foundation=True, + ): + if show: + fig, _axes = self.do_draw(fig_axes, field_max, field_min, title, foundation) + fig.tight_layout() + plt.show() + if save: + _fig, _axes = self.do_draw(fig_axes, field_max, field_min, title, foundation) + self.save_plot(save_format, name=save) + + def do_draw( + self, + fig_axes=None, + field_max=None, + field_min=None, + title=None, + foundation=True, ): fig, axes = fig_axes or plt.subplots() @@ -81,12 +97,7 @@ def draw( axes.set_aspect("equal", adjustable="box") if title is not None: plt.title(title) - - if show: - fig.tight_layout() - plt.show() - if save: - self.save_plot(save_format, name=save) + return fig, axes def set_axes_limits(self, axes, foundation): # pylint: disable=nested-min-max diff --git a/conmech/solvers/optimization/optimization.py b/conmech/solvers/optimization/optimization.py index aa39e122..5ed5b51e 100644 --- a/conmech/solvers/optimization/optimization.py +++ b/conmech/solvers/optimization/optimization.py @@ -205,11 +205,16 @@ def constr(x): sols.append(solution) loss.append(self.loss(solution, *args)[0]) else: + subgrad = None + if method.startswith("gradiented "): + subgrad = self.subgradient + method = method[len("gradiented "):] result = scipy.optimize.minimize( self.loss, solution, args=args, method=method, + jac=subgrad, options={"disp": disp, "maxiter": maxiter}, tol=tol, ) @@ -220,7 +225,6 @@ def constr(x): norm = np.linalg.norm(np.subtract(solution, old_solution)) old_solution = solution.copy() min_index = loss.index(np.min(loss)) - print(method, np.min(loss)) # TODO solution = sols[min_index] return solution diff --git a/examples/Bagirrov_Bartman_Ochal_2024.py b/examples/Bagirrov_Bartman_Ochal_2024.py index cb8378cb..d33906ba 100644 --- a/examples/Bagirrov_Bartman_Ochal_2024.py +++ b/examples/Bagirrov_Bartman_Ochal_2024.py @@ -3,6 +3,7 @@ """ import pickle +import time from dataclasses import dataclass from matplotlib import pyplot as plt @@ -13,7 +14,7 @@ from conmech.properties.mesh_description import RectangleMeshDescription from conmech.scenarios.problems import ContactLaw, StaticDisplacementProblem from conmech.simulations.problem_solver import StaticSolver -from conmech.solvers.optimization.optimization import Optimization +from examples.Makela_et_al_1998 import loss_value, plot_losses mesh_density = 4 kN = 1000 @@ -121,8 +122,11 @@ def main(config: Config, methods, forces): print("Simulating...") setup = StaticSetup(mesh_descr=mesh_descr) + losses = {} for method, force in to_simulate: print(method, force) + m_loss = losses.get(method, {}) + losses[method] = m_loss def outer_forces(x, t=None): if x[1] >= 0.0099: @@ -132,7 +136,9 @@ def outer_forces(x, t=None): setup.outer_forces = outer_forces runner = StaticSolver(setup, "schur") + validator = StaticSolver(setup, "global") + start = time.time() state = runner.solve( verbose=True, fixed_point_abs_tol=0.001, @@ -140,6 +146,7 @@ def outer_forces(x, t=None): method=method, maxiter=100, ) + m_loss[force] = loss_value(state, validator), time.time() - start path = f"{config.outputs_path}/{PREFIX}_mtd_{method}_frc_{force:.2e}" with open(path, "wb+") as output: state.body.dynamics.force.outer.source = None @@ -148,7 +155,15 @@ def outer_forces(x, t=None): state.setup = None state.constitutive_law = None pickle.dump(state, output) - print(Optimization.RESULTS) + path = f"{config.outputs_path}/{PREFIX}_losses" + with open(path, "wb+") as output: + pickle.dump(losses, output) + + print("Plotting...") + + path = f"{config.outputs_path}/{PREFIX}_losses" + if config.show: + plot_losses(path) for m in methods: for f in forces: @@ -185,50 +200,8 @@ def outer_forces(x, t=None): Y[i] = MMLV99().subderivative_normal_direction(X[i], 0.0, 0.0) plt.plot(X, Y) plt.show() - # results = { - # "Powell": [-0.09786211600599237, - # -0.12214289905239312, - # -0.13027212877878766, - # -0.13447218948364842, - # -0.13588717514960513, - # -0.1373096435316275, - # -0.7582249893948801, - # -0.8589012530191608, - # -1.2688709210981735, ], - # "BFGS": [-0.09487765162385353, - # -0.12207326519092926, - # -0.11772218280878324, - # -0.1198269497911567, - # -0.12061095335641955, - # -0.1219350781729528, - # -0.12279409585312624, - # -0.8584230093357013, - # -1.2687124265093166, ], - # "CG": [-0.0955742828809952, - # -0.12191044159984168, - # -0.13009806547436803, - # -0.1341887930175023, - # -0.1358025353476277, - # -0.136904521914724, - # -0.13865495481609302, - # -0.8584104766729636, - # -1.2658836730355307, ], - # "subgradient2": [-0.09786204500205781, - # -0.12214281874382482, - # -0.13027204041320914, - # -0.15450061948841598, - # -0.1571765749815719, - # -0.15986547858657962, - # -0.7582249071621823, - # -0.8589012119331098, - # -1.2688708874747643, ], - # } - # for m, losses in results.items(): - # plt.plot(-1 * np.asarray(losses), label=m) - # plt.legend() - # # plt.loglog() - # plt.show() - methods = ("BFGS", "CG", "qsm", "Powell", "globqsm") + + methods = ("gradiented BFGS", "gradiented CG", "BFGS", "CG", "qsm", "Powell", "globqsm")[:] forces = ( 23e3 * kN, 25e3 * kN, @@ -239,5 +212,5 @@ def outer_forces(x, t=None): 26.2e3 * kN, 27e3 * kN, 30e3 * kN, - )[-1::4] - main(Config(save=False, show=True, force=False).init(), methods, forces) + )[:] + main(Config(save=False, show=False, force=False).init(), methods, forces) diff --git a/examples/Makela_et_al_1998.py b/examples/Makela_et_al_1998.py index fe65f58e..e256ebab 100644 --- a/examples/Makela_et_al_1998.py +++ b/examples/Makela_et_al_1998.py @@ -17,9 +17,11 @@ # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, # USA. import pickle +import time from dataclasses import dataclass import numpy as np +import matplotlib.pyplot as plt from conmech.dynamics.contact.contact_law import PotentialOfContactLaw, DirectContactLaw from conmech.helpers.config import Config @@ -147,8 +149,11 @@ def main(config: Config, methods, forces): print("Simulating...") setup = StaticSetup(mesh_descr=mesh_descr) + losses = {} for method, force in to_simulate: print(method, force) + m_loss = losses.get(method, {}) + losses[method] = m_loss def outer_forces(x, t=None): if x[1] >= 0.0099: @@ -157,14 +162,17 @@ def outer_forces(x, t=None): setup.outer_forces = outer_forces - runner = StaticSolver(setup, "schur") + runner = StaticSolver(setup, "global") + validator = StaticSolver(setup, "global") + start = time.time() state = runner.solve( verbose=True, fixed_point_abs_tol=0.001, initial_displacement=setup.initial_displacement, method=method, ) + m_loss[force] = loss_value(state, validator), time.time() - start path = f"{config.outputs_path}/{PREFIX}_mtd_{method}_frc_{force:.2e}" with open(path, "wb+") as output: state.body.dynamics.force.outer.source = None @@ -173,7 +181,15 @@ def outer_forces(x, t=None): state.setup = None state.constitutive_law = None pickle.dump(state, output) - print(Optimization.RESULTS) + path = f"{config.outputs_path}/{PREFIX}_losses" + with open(path, "wb+") as output: + pickle.dump(losses, output) + + print("Plotting...")\ + + path = f"{config.outputs_path}/{PREFIX}_losses" + if config.show: + plot_losses(path) for m in methods: for f in forces: @@ -181,7 +197,6 @@ def outer_forces(x, t=None): with open(path, "rb") as output: state = pickle.load(output) - print("drawing") drawer = Drawer(state=state, config=config) drawer.colorful = True drawer.draw( @@ -190,6 +205,7 @@ def outer_forces(x, t=None): # title=f"{m}: {f}, " # f"time: {runner.step_solver.last_timing}" ) + x = state.body.mesh.nodes[: state.body.mesh.contact_nodes_count - 1, 0] u = state.displacement[: state.body.mesh.contact_nodes_count - 1, 1] y1 = [MMLV99().subderivative_normal_direction(-u_, 0.0, 0.0) for u_ in u] @@ -202,9 +218,53 @@ def outer_forces(x, t=None): plt.show() -if __name__ == "__main__": - from matplotlib import pyplot as plt +def plot_losses(path): + print("Plotting...") + with open(path, "rb") as output: + losses = pickle.load(output) + + for mtd, values in losses.items(): + forces_ = np.asarray(list(values.keys())) + values_ = np.asarray(list(values.values()))[:, 0] + times_ = np.asarray(list(values.values()))[:, 1] + plt.plot(forces_ / 1e3, -1 * values_, "-o", label=mtd) + plt.legend() + plt.ylabel("$-\mathcal{L}(u)$") + plt.xlabel(r"Load [kN/m$^2$]") + plt.grid() + plt.show() + for mtd, values in losses.items(): + forces_ = np.asarray(list(values.keys())) + values_ = np.asarray(list(values.values()))[:, 0] + times_ = np.asarray(list(values.values()))[:, 1] + plt.plot(forces_ / 1e3, times_, "-o", label=mtd) + plt.legend() + plt.ylabel("$time$") + plt.xlabel(r"Load [kN/m$^2$]") + plt.grid() + plt.show() + +def loss_value(state, runner) -> float: + initial_guess = state["displacement"].T.ravel().reshape(state.body.mesh.dimension, -1) + solution = np.squeeze(initial_guess.copy().reshape(1, -1)) + self: Optimization = runner.step_solver + args = ( + np.zeros_like(solution), # variable_old + self.body.mesh.nodes, + self.body.mesh.contact_boundary, + self.body.mesh.boundaries.contact_normals, + self.lhs, + self.rhs, + np.zeros_like(solution), # displacement + np.ascontiguousarray(self.body.dynamics.acceleration_operator.SM1.data), + self.time_step, + ) + result = self.loss(solution, *args)[0] + print(result) + return result + +if __name__ == "__main__": # X = np.linspace(0, -3 * mm, 1000) # Y = np.empty(1000) # for i in range(1000): @@ -215,44 +275,7 @@ def outer_forces(x, t=None): # Y[i] = MMLV99.normal_direction(X[i]) # plt.plot(X, Y) # plt.show() - # results = { - # "BFGS": [-0.061546678021737036, - # -0.06782602334922566, - # -0.07441012406759984, - # -0.08129875924227234, - # -0.0959892642846613, - # -0.10379118250601398, - # -0.10538811134540409, - # -0.8584224789292736, - # -0.14133884664811114, ], - # "CG": [-0.07225702623584927, - # -0.07966800277816762, - # -0.08744039267159345, - # -0.09557428287965247, - # -0.12191044159984168, - # -0.1358025353476277, - # -0.13865495481609302, - # -0.15028696247286885, - # -1.265832916470563, ], - # "Powell": [-0.0723012449592487, - # -0.07971212256709342, - # -0.0874845064006726, - # -0.0978621160055679, - # -0.12214289905071576, - # -0.13588717513833654, - # -0.7582249892835198, - # -0.8589012526317955, - # -1.2688709207679356, ], - # "subgradient": [-0.05079652409797247, - # -0.046161334145372934, - # -0.04120648554585715, - # -0.3859157295854724, - # -0.6104716467978587, - # -0.7302821710666211, - # -0.7554950402698594, - # -0.8555741662642888, - # -1.2663638426265278, ], - # } + forces = np.asarray( ( 20e3 * kN, @@ -263,15 +286,9 @@ def outer_forces(x, t=None): 26e3 * kN, 26.2e3 * kN, 27e3 * kN, - 30e3 * kN, + # 30e3 * kN, ) - )[::2] - # # for m, losses in results.items(): - # # plt.plot(forces/1e3, -1 * np.asarray(losses[:]), "-o", label=m) - # plt.legend() - # plt.ylabel("$-\mathcal{L}(u)$") - # plt.xlabel(r"Load [kN/m$^2$]") - # plt.grid() - # plt.show() - methods = ("BFGS", "CG", "Powell", "globqsm")[-2:] - main(Config(save=False, show=True, force=True).init(), methods, forces) + )[:] + + methods = ("gradiented BFGS", "gradiented CG", "BFGS", "CG", "Powell", "globqsm", "qsm")[:] + main(Config(save=False, show=False, force=True).init(), methods, forces) From e6f3b7f11b481137cdb3890db367fed718a1093c Mon Sep 17 00:00:00 2001 From: Piotr Bartman-Szwarc Date: Thu, 31 Oct 2024 10:33:25 +0100 Subject: [PATCH 29/34] BBO: DC functions --- conmech/solvers/optimization/optimization.py | 9 +++- conmech/solvers/solver_methods.py | 6 ++- examples/Bagirrov_Bartman_Ochal_2024.py | 25 +++++++++++ examples/Makela_et_al_1998.py | 46 +++++++++++++++++--- 4 files changed, 77 insertions(+), 9 deletions(-) diff --git a/conmech/solvers/optimization/optimization.py b/conmech/solvers/optimization/optimization.py index 5ed5b51e..7e7e548b 100644 --- a/conmech/solvers/optimization/optimization.py +++ b/conmech/solvers/optimization/optimization.py @@ -84,6 +84,13 @@ def __init__( ) else: self.subgradient = None + if hasattr(contact_law, "sub2derivative_normal_direction"): # TODO + self.sub2gradient = make_subgradient( + normal_condition=contact_law.sub2derivative_normal_direction, + only_boundary=True, + ) + else: + self.sub2gradient = None if isinstance(statement, WaveStatement): if isinstance(contact_law, InteriorContactLaw): self.loss = make_equation( # TODO! @@ -150,7 +157,7 @@ def _solve_impl( # pylint: disable=import-outside-toplevel,import-error) from kosopt.qsmlm import make_minimizer - self.minimizer = make_minimizer(self.loss, self.subgradient) + self.minimizer = make_minimizer(self.loss, self.subgradient, self.sub2gradient) while norm >= fixed_point_abs_tol: if method.lower() in QSMLM_NAMES: diff --git a/conmech/solvers/solver_methods.py b/conmech/solvers/solver_methods.py index dfb2e256..df857d90 100644 --- a/conmech/solvers/solver_methods.py +++ b/conmech/solvers/solver_methods.py @@ -228,6 +228,7 @@ def make_subgradient( tangential_condition_bound: Optional[Callable] = None, problem_dimension=2, variable_dimension=2, + only_boundary=False, ): normal_condition = njit(normal_condition) normal_condition_bound = njit(normal_condition_bound, value=1) @@ -304,7 +305,10 @@ def cost_functional( var, var_old, u_vector, nodes, contact_boundary, contact_normals, dt ) ind = lhs.shape[0] - result_ = np.dot(lhs, var[:ind]) - rhs + dj + if only_boundary: + result_ = dj + else: + result_ = np.dot(lhs, var[:ind]) - rhs + dj result[:ind] = result_.ravel() return result diff --git a/examples/Bagirrov_Bartman_Ochal_2024.py b/examples/Bagirrov_Bartman_Ochal_2024.py index d33906ba..0513ba2a 100644 --- a/examples/Bagirrov_Bartman_Ochal_2024.py +++ b/examples/Bagirrov_Bartman_Ochal_2024.py @@ -65,6 +65,27 @@ def subderivative_normal_direction( return k20 * (var_nu * 2) + k21 return var_nu**3 * 4 * k30 + @staticmethod + def sub2derivative_normal_direction( + var_nu: float, static_displacement_nu: float, dt: float + ) -> float: + acc = 0.0 + p1 = 0.5 * mm + if var_nu <= p1: + return acc + + acc += (k0 * p1 * 2) - (k10 * (p1 * 2) + k11) + p2 = 1 * mm + if var_nu < p2: + return acc + + acc += (k10 * (p2 * 2) + k11) - (k20 * (p2 * 2) + k21) + p3 = 2 * mm + if var_nu < p3: + return acc + + return acc + @dataclass() class StaticSetup(StaticDisplacementProblem): @@ -200,6 +221,10 @@ def outer_forces(x, t=None): Y[i] = MMLV99().subderivative_normal_direction(X[i], 0.0, 0.0) plt.plot(X, Y) plt.show() + for i in range(1000): + Y[i] = MMLV99().sub2derivative_normal_direction(X[i], 0.0, 0.0) + plt.plot(X, Y) + plt.show() methods = ("gradiented BFGS", "gradiented CG", "BFGS", "CG", "qsm", "Powell", "globqsm")[:] forces = ( diff --git a/examples/Makela_et_al_1998.py b/examples/Makela_et_al_1998.py index e256ebab..dd3841d5 100644 --- a/examples/Makela_et_al_1998.py +++ b/examples/Makela_et_al_1998.py @@ -76,6 +76,30 @@ def subderivative_normal_direction( return k20 * (u_nu * 2) + k21 return 0.0 + @staticmethod + def sub2derivative_normal_direction( + var_nu: float, static_displacement_nu: float, dt: float + ) -> float: + u_nu = -var_nu + + acc = 0.0 + p1 = 0.5 * mm + if u_nu <= p1: + return acc + + acc += (k0 * p1 * 2) - (k10 * (p1 * 2) + k11) + p2 = 1.0 * mm + if u_nu < p2: + return acc + + acc += (k10 * (p2 * 2) + k11) - (k20 * (p2 * 2) + k21) + p3 = 2.0 * mm + if u_nu < p3: + return acc + + acc += (k20 * (p3 * 2) + k21) - 0.0 + return acc + @staticmethod def potential_tangential_direction( var_tau: float, static_displacement_tau: float, dt: float @@ -188,8 +212,8 @@ def outer_forces(x, t=None): print("Plotting...")\ path = f"{config.outputs_path}/{PREFIX}_losses" - if config.show: - plot_losses(path) + # if config.show: + plot_losses(path) for m in methods: for f in forces: @@ -265,14 +289,22 @@ def loss_value(state, runner) -> float: if __name__ == "__main__": - # X = np.linspace(0, -3 * mm, 1000) - # Y = np.empty(1000) + X = np.linspace(0, -3 * mm, 1000) + Y = np.empty(1000) + for i in range(1000): + Y[i] = MMLV99.potential_normal_direction(X[i], 0, 0) + plt.plot(X, Y) + plt.show() + for i in range(1000): + Y[i] = MMLV99.subderivative_normal_direction(X[i], 0, 0) + plt.plot(X, Y) + plt.show() # for i in range(1000): - # Y[i] = MMLV99.potential_normal_direction(X[i]) + # Y[i] = MMLV99.sub2derivative_normal_direction(X[i], 0, 0) # plt.plot(X, Y) # plt.show() # for i in range(1000): - # Y[i] = MMLV99.normal_direction(X[i]) + # Y[i] = MMLV99.subderivative_normal_direction(X[i], 0, 0) + MMLV99.sub2derivative_normal_direction(X[i], 0, 0) # plt.plot(X, Y) # plt.show() @@ -290,5 +322,5 @@ def loss_value(state, runner) -> float: ) )[:] - methods = ("gradiented BFGS", "gradiented CG", "BFGS", "CG", "Powell", "globqsm", "qsm")[:] + methods = ("gradiented BFGS", "gradiented CG", "BFGS", "CG", "Powell", "globqsm", "qsm")[-1:] main(Config(save=False, show=False, force=True).init(), methods, forces) From 4d46d3eab7099f36a0cf1a5d1347044bf03268c7 Mon Sep 17 00:00:00 2001 From: Piotr Bartman-Szwarc Date: Tue, 5 Nov 2024 17:24:19 +0100 Subject: [PATCH 30/34] BBO: dc as option --- conmech/solvers/optimization/optimization.py | 8 ++++++- examples/Makela_et_al_1998.py | 24 ++++++++++++-------- 2 files changed, 21 insertions(+), 11 deletions(-) diff --git a/conmech/solvers/optimization/optimization.py b/conmech/solvers/optimization/optimization.py index 7e7e548b..82b7b1e8 100644 --- a/conmech/solvers/optimization/optimization.py +++ b/conmech/solvers/optimization/optimization.py @@ -45,6 +45,7 @@ "qsm", "qsmlm", } +QSMLM_NAMES = {"dc " + name for name in QSMLM_NAMES}.union(QSMLM_NAMES) GLOBAL_QSMLM_NAMES = { "global quasi secant method", "global limited memory quasi secant method", @@ -52,6 +53,7 @@ "globqsm", "globqsmlm", } +GLOBAL_QSMLM_NAMES = {"dc " + name for name in GLOBAL_QSMLM_NAMES}.union(GLOBAL_QSMLM_NAMES) class Optimization(Solver): @@ -157,7 +159,11 @@ def _solve_impl( # pylint: disable=import-outside-toplevel,import-error) from kosopt.qsmlm import make_minimizer - self.minimizer = make_minimizer(self.loss, self.subgradient, self.sub2gradient) + self.minimizer = make_minimizer( + self.loss, + self.subgradient, + self.sub2gradient if method.lower().startswith("dc") else None + ) while norm >= fixed_point_abs_tol: if method.lower() in QSMLM_NAMES: diff --git a/examples/Makela_et_al_1998.py b/examples/Makela_et_al_1998.py index dd3841d5..4f0806e3 100644 --- a/examples/Makela_et_al_1998.py +++ b/examples/Makela_et_al_1998.py @@ -247,11 +247,15 @@ def plot_losses(path): with open(path, "rb") as output: losses = pickle.load(output) - for mtd, values in losses.items(): + for i, (mtd, values) in enumerate(losses.items()): forces_ = np.asarray(list(values.keys())) values_ = np.asarray(list(values.values()))[:, 0] times_ = np.asarray(list(values.values()))[:, 1] - plt.plot(forces_ / 1e3, -1 * values_, "-o", label=mtd) + # plt.plot(forces_ / 1e3, -1 * values_, "-o", label=mtd, linewidth=3*(len(losses.items()) - i)) + total_width = 300 + width = total_width / len(losses.items()) + shift = total_width / 2 + i * width + width / 2 + plt.bar(forces_ / 1e3 + shift, -1 * values_, label=mtd, width=width) plt.legend() plt.ylabel("$-\mathcal{L}(u)$") plt.xlabel(r"Load [kN/m$^2$]") @@ -299,14 +303,14 @@ def loss_value(state, runner) -> float: Y[i] = MMLV99.subderivative_normal_direction(X[i], 0, 0) plt.plot(X, Y) plt.show() - # for i in range(1000): - # Y[i] = MMLV99.sub2derivative_normal_direction(X[i], 0, 0) - # plt.plot(X, Y) - # plt.show() - # for i in range(1000): - # Y[i] = MMLV99.subderivative_normal_direction(X[i], 0, 0) + MMLV99.sub2derivative_normal_direction(X[i], 0, 0) - # plt.plot(X, Y) - # plt.show() + for i in range(1000): + Y[i] = MMLV99.sub2derivative_normal_direction(X[i], 0, 0) + plt.plot(X, Y) + plt.show() + for i in range(1000): + Y[i] = MMLV99.subderivative_normal_direction(X[i], 0, 0) + MMLV99.sub2derivative_normal_direction(X[i], 0, 0) + plt.plot(X, Y) + plt.show() forces = np.asarray( ( From 51bb4d612d17d16253d7c1f21ec6663f379ae30c Mon Sep 17 00:00:00 2001 From: Piotr Bartman-Szwarc Date: Tue, 5 Nov 2024 17:24:55 +0100 Subject: [PATCH 31/34] BBO:WIP: example 1 --- examples/Bagirrov_Bartman_Ochal_2024.py | 198 +++++++++++++++++++----- 1 file changed, 158 insertions(+), 40 deletions(-) diff --git a/examples/Bagirrov_Bartman_Ochal_2024.py b/examples/Bagirrov_Bartman_Ochal_2024.py index 0513ba2a..95220596 100644 --- a/examples/Bagirrov_Bartman_Ochal_2024.py +++ b/examples/Bagirrov_Bartman_Ochal_2024.py @@ -87,6 +87,73 @@ def sub2derivative_normal_direction( return acc +def make_composite_contact_law(layers, alpha, beta): + alphas = np.asarray([alpha(lay) for lay in layers]) + betas = np.asarray([beta(lay) for lay in layers]) + betas[0] = 0 + integ = np.asarray([0.5 * (layers[n+1] - layers[n]) * (betas[n+1] + alphas[n]) + for n in range(len(layers) - 1)]) + cumm = np.cumsum(integ) + slope = np.asarray([(betas[n + 1] - alphas[n]) / (layers[n + 1] - layers[n]) for n in range(len(layers) - 1)]) + const = np.asarray([alphas[n] - layers[n] * slope[n] for n in range(len(layers) - 1)]) + + class Composite(ContactLaw): + + @staticmethod + def potential_tangential_direction( + var_tau: float, static_displacement_tau: float, dt: float + ) -> float: + return np.log(np.sum(var_tau ** 2) ** 0.5 + 1) + + @staticmethod + def potential_normal_direction( + var_nu: float, static_displacement_nu: float, dt: float + ) -> float: + val = 0.0 + if var_nu <= layers[0]: + return 0.0 + for n in range(len(layers)): + if var_nu > layers[n]: + continue + if n > 1: + val += cumm[n-2] + val += 0.5 * (var_nu - layers[n-1]) * (var_nu * slope[n-1] + const[n-1] + alphas[n-1]) + break + else: + val = cumm[-1] + alphas[-1] * (var_nu - layers[-1]) + return val + + @staticmethod + def subderivative_normal_direction( + var_nu: float, static_displacement_nu: float, dt: float + ) -> float: + if var_nu <= layers[0]: + return 0.0 + if layers[-1] < var_nu: + return alphas[-1] + for n in range(len(layers) - 1): + if var_nu <= layers[n+1]: + return var_nu * slope[n] + const[n] + + @staticmethod + def sub2derivative_normal_direction( + var_nu: float, static_displacement_nu: float, dt: float + ) -> float: + acc = 0.0 + if var_nu <= layers[0]: + return 0.0 + for n in range(len(layers)): + if var_nu >= layers[n]: + acc += max(0, betas[n] - alphas[n]) + if n > 0: + acc += max(0, alphas[n-1] - betas[n]) + else: + acc += max(0, - slope[n-1]) * (var_nu - layers[n-1]) + break + return acc + + return Composite + @dataclass() class StaticSetup(StaticDisplacementProblem): grid_height: ... = 10 * mm @@ -112,7 +179,7 @@ def friction_bound(u_nu: float) -> float: ) -def main(config: Config, methods, forces): +def main(config: Config, methods, forces, contact): """ Entrypoint to example. @@ -141,7 +208,7 @@ def main(config: Config, methods, forces): if to_simulate: print("Simulating...") - setup = StaticSetup(mesh_descr=mesh_descr) + setup = StaticSetup(mesh_descr=mesh_descr, contact_law=contact) losses = {} for method, force in to_simulate: @@ -183,59 +250,110 @@ def outer_forces(x, t=None): print("Plotting...") path = f"{config.outputs_path}/{PREFIX}_losses" - if config.show: - plot_losses(path) - - for m in methods: - for f in forces: - path = f"{config.outputs_path}/{PREFIX}_mtd_{m}_frc_{f:.2e}" - with open(path, "rb") as output: - state = pickle.load(output) - - drawer = Drawer(state=state, config=config) - drawer.colorful = True - drawer.draw( - show=config.show, - save=config.save, - title=f"{m}: {f}, ", - # f"time: {runner.step_solver.last_timing}" - ) - x = state.body.mesh.nodes[: state.body.mesh.contact_nodes_count - 1, 0] - u = state.displacement[: state.body.mesh.contact_nodes_count - 1, 1] - y1 = [MMLV99().subderivative_normal_direction(-u_, 0.0, 0.0) for u_ in u] - print(f) - plt.plot(x, y1, label=f"{f:.2e}") - plt.title(m) - plt.legend() - plt.show() + # if config.show: + plot_losses(path) + + # for m in methods: + # for f in forces: + # path = f"{config.outputs_path}/{PREFIX}_mtd_{m}_frc_{f:.2e}" + # with open(path, "rb") as output: + # state = pickle.load(output) + # + # drawer = Drawer(state=state, config=config) + # drawer.colorful = True + # drawer.draw( + # show=config.show, + # save=config.save, + # title=f"{m}: {f}, ", + # # f"time: {runner.step_solver.last_timing}" + # ) + # x = state.body.mesh.nodes[: state.body.mesh.contact_nodes_count - 1, 0] + # u = state.displacement[: state.body.mesh.contact_nodes_count - 1, 1] + # y1 = [MMLV99().subderivative_normal_direction(-u_, 0.0, 0.0) for u_ in u] + # print(f) + # plt.plot(x, y1, label=f"{f:.2e}") + # plt.title(m) + # plt.legend() + # plt.show() + + +def composite_problem(config, layers_limits, thickness, methods, forces): + def bottom(x): + if x >= thickness: + return 0.0 + return 0.0 #(x / mm) ** 2 / mm + def top(x): + # if x >= thickness: + # return 0.0 + return 15.0 / mm + (x / mm) / mm -if __name__ == "__main__": - X = np.linspace((2 - 2) * mm, (2 + 2) * mm, 1000) + contact = make_composite_contact_law(layers_limits, alpha=bottom, beta=top) + + X = np.linspace(-thickness / 4, 5 / 4 * thickness, 1000) Y = np.empty(1000) for i in range(1000): - Y[i] = MMLV99.potential_normal_direction(X[i], 0.0, 0.0) + Y[i] = contact.potential_normal_direction(X[i], 0.0, 0.0) plt.plot(X, Y) plt.show() for i in range(1000): - Y[i] = MMLV99().subderivative_normal_direction(X[i], 0.0, 0.0) + Y[i] = contact.subderivative_normal_direction(X[i], 0.0, 0.0) plt.plot(X, Y) - plt.show() + # plt.show() for i in range(1000): - Y[i] = MMLV99().sub2derivative_normal_direction(X[i], 0.0, 0.0) + Y[i] = contact.sub2derivative_normal_direction(X[i], 0.0, 0.0) plt.plot(X, Y) plt.show() - - methods = ("gradiented BFGS", "gradiented CG", "BFGS", "CG", "qsm", "Powell", "globqsm")[:] + main(config, methods, forces, contact) + + +def survey(config): + methods = ( + "gradiented BFGS", + # "gradiented CG", + "BFGS", + # "CG", + "Powell", + "qsm", + "dc qsm", + "globqsm", + "dc globqsm" + )[:] forces = ( 23e3 * kN, 25e3 * kN, - 25.6e3 * kN, - 25.9e3 * kN, + 25.5e3 * kN, 26e3 * kN, - 26.1e3 * kN, - 26.2e3 * kN, + 26.5e3 * kN, 27e3 * kN, + 28e3 * kN, 30e3 * kN, )[:] - main(Config(save=False, show=False, force=False).init(), methods, forces) + for i in [6] : #range(1, 10 + 1, 2): + thickness = 3 * mm + layers_num = i + # layers_limits = np.logspace(0, thickness, layers_num + 1) + layers_limits = partition(0, thickness, layers_num + 1, p=1.25) + composite_problem(config, layers_limits, thickness, methods, forces) + + +def partition(start, stop, num, p=1.0): + length = stop - start + if p == 1: + grad = num + else: + grad = (1 - p**num) / (1 - p) + + first_part = length / grad + points = [start] + curr = start + for i in range(num): + curr += first_part * (p ** i) + points.append(curr) + + return np.asarray(points) + + +if __name__ == "__main__": + config_ = Config(save=False, show=False, force=True).init() + survey(config_) From 31f9e9d97309471263cb589244ecae1bc9f3d502 Mon Sep 17 00:00:00 2001 From: Piotr Bartman-Szwarc Date: Wed, 6 Nov 2024 12:41:13 +0100 Subject: [PATCH 32/34] BBO:WIP: example 1 --- examples/Bagirrov_Bartman_Ochal_2024.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/Bagirrov_Bartman_Ochal_2024.py b/examples/Bagirrov_Bartman_Ochal_2024.py index 95220596..2c506115 100644 --- a/examples/Bagirrov_Bartman_Ochal_2024.py +++ b/examples/Bagirrov_Bartman_Ochal_2024.py @@ -179,13 +179,13 @@ def friction_bound(u_nu: float) -> float: ) -def main(config: Config, methods, forces, contact): +def main(config: Config, methods, forces, contact, prefix=""): """ Entrypoint to example. To see result of simulation you need to call from python `main(Config().init())`. """ - PREFIX = "DBBO" + PREFIX = f"BBO{prefix}" if config.force: to_simulate = [(m, f) for m in methods for f in forces] else: @@ -304,7 +304,7 @@ def top(x): Y[i] = contact.sub2derivative_normal_direction(X[i], 0.0, 0.0) plt.plot(X, Y) plt.show() - main(config, methods, forces, contact) + main(config, methods, forces, contact, prefix=str(len(layers_limits))) def survey(config): @@ -329,7 +329,7 @@ def survey(config): 28e3 * kN, 30e3 * kN, )[:] - for i in [6] : #range(1, 10 + 1, 2): + for i in [0, 1, 2, 3, 4, 5, 6, 7, 8]: #range(1, 10 + 1, 2): thickness = 3 * mm layers_num = i # layers_limits = np.logspace(0, thickness, layers_num + 1) From 2d80adab830fd534e7fe4f5902ca5efef78f05b6 Mon Sep 17 00:00:00 2001 From: Piotr Bartman-Szwarc Date: Tue, 12 Nov 2024 12:17:38 +0100 Subject: [PATCH 33/34] BBO:WIP: example 2 --- conmech/scenarios/problems.py | 2 +- conmech/simulations/problem_solver.py | 2 + conmech/solvers/optimization/optimization.py | 28 +++- conmech/solvers/solver.py | 5 +- conmech/solvers/solver_methods.py | 119 +++++++++++++- conmech/struct/types.py | 30 ++++ examples/Bagirrov_Bartman_Ochal_2024.py | 157 ++++++------------- examples/Makela_et_al_1998.py | 4 +- 8 files changed, 220 insertions(+), 127 deletions(-) create mode 100644 conmech/struct/types.py diff --git a/conmech/scenarios/problems.py b/conmech/scenarios/problems.py index c9589705..bd996fb5 100644 --- a/conmech/scenarios/problems.py +++ b/conmech/scenarios/problems.py @@ -62,7 +62,7 @@ class StaticProblem(Problem, ABC): class TimeDependentProblem(Problem, ABC): - time_step = 0 + time_step = 0.0 class QuasistaticProblem(TimeDependentProblem, ABC): diff --git a/conmech/simulations/problem_solver.py b/conmech/simulations/problem_solver.py index 13d9037d..190dca39 100644 --- a/conmech/simulations/problem_solver.py +++ b/conmech/simulations/problem_solver.py @@ -84,6 +84,8 @@ def __init__(self, problem: Problem, body_properties: BodyProperties): self.done = 0 self.to_do = 1 + self.computation_time = None + @property def solving_method(self) -> str: return str(self.step_solver) diff --git a/conmech/solvers/optimization/optimization.py b/conmech/solvers/optimization/optimization.py index 82b7b1e8..211ead97 100644 --- a/conmech/solvers/optimization/optimization.py +++ b/conmech/solvers/optimization/optimization.py @@ -17,12 +17,14 @@ # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, # USA. import math +import time from typing import Optional import numba import numpy as np import scipy.optimize +from conmech.struct.types import * from conmech.dynamics.statement import ( Statement, StaticPoissonStatement, @@ -34,7 +36,7 @@ from conmech.solvers.solver_methods import ( make_cost_functional, make_equation, - make_subgradient, + make_subgradient, make_subgradient_dc, ) @@ -87,9 +89,10 @@ def __init__( else: self.subgradient = None if hasattr(contact_law, "sub2derivative_normal_direction"): # TODO - self.sub2gradient = make_subgradient( - normal_condition=contact_law.sub2derivative_normal_direction, - only_boundary=True, + self.sub2gradient = make_subgradient_dc( + normal_condition=contact_law.subderivative_normal_direction, + normal_condition_sub2=contact_law.sub2derivative_normal_direction, + # only_boundary=True, ) else: self.sub2gradient = None @@ -131,6 +134,7 @@ def _solve_impl( fixed_point_abs_tol: float = math.inf, **kwargs, ) -> np.ndarray: + start = time.time() norm = math.inf solution = np.squeeze(initial_guess.copy().reshape(1, -1)) displacement = np.squeeze(displacement.copy().reshape(1, -1)) @@ -144,7 +148,7 @@ def _solve_impl( self.body.mesh.contact_boundary, self.body.mesh.boundaries.contact_normals, self.lhs, - self.rhs, + self.rhs if len(self.rhs.shape) == 1 else self.rhs[0], # TODO displacement, np.ascontiguousarray(self.body.dynamics.acceleration_operator.SM1.data), self.time_step, @@ -154,6 +158,7 @@ def _solve_impl( sols = [] sols.append(solution) loss.append(self.loss(solution, *args)[0]) + self.computation_time += time.time() - start if self.minimizer is None and method.lower() in QSMLM_NAMES.union(GLOBAL_QSMLM_NAMES): # pylint: disable=import-outside-toplevel,import-error) @@ -167,23 +172,27 @@ def _solve_impl( while norm >= fixed_point_abs_tol: if method.lower() in QSMLM_NAMES: - solution = self.minimizer(solution, args, maxiter=maxiter) + start = time.time() + solution = self.minimizer(solution, args, 0, 1, maxiter) + self.computation_time += time.time() - start sols.append(solution.copy()) loss.append(self.loss(solution, *args)[0]) elif method.lower() in GLOBAL_QSMLM_NAMES: # pylint: disable=import-outside-toplevel,import-error) from kosopt import subgradient - solution = subgradient.minimize( + solution, comp_time = subgradient.minimize( self.minimizer, self.loss, solution, args, maxiter=maxiter, subgradient=self.subgradient, + sub2gradient=self.sub2gradient, ) sols.append(solution.copy()) loss.append(self.loss(solution, *args)[0]) + self.computation_time += comp_time elif method.lower() in ( # TODO "discontinuous gradient", "discontinuous gradient method", @@ -192,7 +201,7 @@ def _solve_impl( # pylint: disable=import-outside-toplevel,import-error) from kosopt import qsmlmi - solution = qsmlmi.minimize(self.loss, solution, args=args, maxiter=maxiter) + solution = qsmlmi.minimize(self.loss, solution, args, 0, 1, 10000) sols.append(solution.copy()) loss.append(self.loss(solution, *args)[0]) elif method.lower() == "constrained": @@ -222,6 +231,7 @@ def constr(x): if method.startswith("gradiented "): subgrad = self.subgradient method = method[len("gradiented "):] + start = time.time() result = scipy.optimize.minimize( self.loss, solution, @@ -231,12 +241,14 @@ def constr(x): options={"disp": disp, "maxiter": maxiter}, tol=tol, ) + self.computation_time += time.time() - start solution = result.x sols.append(solution.copy()) loss.append(self.loss(solution, *args)[0]) norm = np.linalg.norm(np.subtract(solution, old_solution)) old_solution = solution.copy() + break min_index = loss.index(np.min(loss)) solution = sols[min_index] diff --git a/conmech/solvers/solver.py b/conmech/solvers/solver.py index 296278c4..378019b0 100644 --- a/conmech/solvers/solver.py +++ b/conmech/solvers/solver.py @@ -63,7 +63,7 @@ def __init__( ) ) - self.last_timing = None + self.computation_time = 0.0 def __str__(self) -> str: raise NotImplementedError() @@ -82,7 +82,6 @@ def _solve_impl( raise NotImplementedError() def solve(self, initial_guess: np.ndarray, **kwargs) -> np.ndarray: - start = time.time() solution = self._solve_impl( initial_guess, variable_old=self.v_vector, @@ -98,6 +97,4 @@ def solve(self, initial_guess: np.ndarray, **kwargs) -> np.ndarray: ): solution[i] = c[j] - self.last_timing = time.time() - start - return solution diff --git a/conmech/solvers/solver_methods.py b/conmech/solvers/solver_methods.py index df857d90..5d329b96 100644 --- a/conmech/solvers/solver_methods.py +++ b/conmech/solvers/solver_methods.py @@ -21,9 +21,11 @@ import numba import numpy as np +from conmech.struct.types import * from conmech.helpers import nph + @numba.njit(inline="always") def interpolate_node_between(edge, vector, full_vector, dimension): result = np.zeros(dimension) @@ -153,11 +155,11 @@ def make_cost_functional( tangential_condition = njit(tangential_condition) tangential_condition_bound = njit(tangential_condition_bound, value=1) - @numba.njit() + @numba.njit(f64(f64, f64, f64, f64, f64)) def contact_cost(length, normal, normal_bound, tangential, tangential_bound): return length * (normal_bound * normal + tangential_bound * tangential) - @numba.njit() + @numba.njit(f64(f64[:], const_f64_vec, const_f64_vec, const_f64_mat, const_i64_mat, const_f64_mat, f64)) def contact_cost_functional( var, var_old, static_displacement, nodes, contact_boundary, contact_normals, dt ): @@ -197,7 +199,9 @@ def contact_cost_functional( return cost # pylint: disable=too-many-arguments,unused-argument # 'base_integrals' - @numba.njit() + # @numba.njit(f64[:](f64[:], f64[:], f64[:,:], i64[:,:], f64[:,:], f64[:,:], f64[:], f64[:], f64[:,:], i64)) + @numba.njit(f64[:](f64[:], const_f64_vec, const_f64_mat, const_i64_mat, const_f64_mat, const_f64_mat, const_f64_vec, const_f64_vec, const_f64_mat, i64)) + # @numba.njit() def cost_functional( var, var_old, @@ -313,3 +317,112 @@ def cost_functional( return result return cost_functional + + +def make_subgradient_dc( + normal_condition: Callable, + normal_condition_sub2: Callable, + normal_condition_bound: Optional[Callable] = None, + tangential_condition: Optional[Callable] = None, + tangential_condition_bound: Optional[Callable] = None, + problem_dimension=2, + variable_dimension=2, + only_boundary=False, +): + normal_condition = njit(normal_condition) + normal_condition_sub2 = njit(normal_condition_sub2) + normal_condition_bound = njit(normal_condition_bound, value=1) + tangential_condition = njit(tangential_condition) + tangential_condition_bound = njit(tangential_condition_bound, value=1) + + @numba.njit() + def contact_cost(length, normal, normal_bound, tangential, tangential_bound): + return length * (normal_bound * normal + tangential_bound * tangential) + + @numba.njit() + def contact_subgradient( + var, var1, var_old, static_displacement, nodes, contact_boundary, contact_normals, dt + ): + cost = np.zeros_like(var) + offset = len(var) // variable_dimension + + # pylint: disable=not-an-iterable + for ei in numba.prange(len(contact_boundary)): + edge = contact_boundary[ei] + normal_vector = contact_normals[ei] + # ASSUMING `u_vector` and `nodes` have the same order! + vm = interpolate_node_between(edge, var, var_old, dimension=variable_dimension) + vm1 = interpolate_node_between(edge, var1, var_old, dimension=variable_dimension) + if variable_dimension == 1: + raise NotImplementedError() # TODO + # vm_normal = vm[0] + # vm_tangential = np.empty(0) + # else: + vm_normal = (vm * normal_vector).sum() + vm_tangential = vm - vm_normal * normal_vector + vm_normal1 = (vm1 * normal_vector).sum() + vm_tangential1 = vm1 - vm_normal1 * normal_vector + + static_displacement_mean = interpolate_node_between( + edge, + static_displacement, + static_displacement, + dimension=problem_dimension, + ) + static_displacement_normal = (static_displacement_mean * normal_vector).sum() + static_displacement_tangential = ( + static_displacement_mean - static_displacement_normal * normal_vector + ) + + subgrad = contact_cost( + nph.length(edge, nodes), + normal_condition(vm_normal1, static_displacement_normal, dt), + normal_condition_bound(vm_normal1, static_displacement_normal, dt), + tangential_condition(vm_tangential, static_displacement_tangential, dt), + tangential_condition_bound(vm_normal1, static_displacement_normal, dt), + ) + contact_cost( + nph.length(edge, nodes), + normal_condition_sub2(vm_normal1, static_displacement_normal, dt), + normal_condition_bound(vm_normal1, static_displacement_normal, dt), + tangential_condition(vm_tangential, static_displacement_tangential, dt), + tangential_condition_bound(vm_normal1, static_displacement_normal, dt), + ) - contact_cost( + nph.length(edge, nodes), + normal_condition_sub2(vm_normal, static_displacement_normal, dt), + normal_condition_bound(vm_normal, static_displacement_normal, dt), + tangential_condition(vm_tangential, static_displacement_tangential, dt), + tangential_condition_bound(vm_normal, static_displacement_normal, dt), + ) + + for node in edge: + for i in range(variable_dimension): + if node < offset: + cost[i * offset + node] += normal_vector[i] / len(edge) * subgrad + + return cost + + # pylint: disable=too-many-arguments,unused-argument # 'base_integrals' + @numba.njit() + def cost_functional( + var, + var1, + var_old, + nodes, + contact_boundary, + contact_normals, + lhs, + rhs, + u_vector, + base_integrals, + dt, + ): + result = np.zeros_like(var) + dj = contact_subgradient( + var, var1, var_old, u_vector, nodes, contact_boundary, contact_normals, dt + ) + ind = lhs.shape[0] + result_ = np.dot(lhs, var1[:ind]) - rhs + dj + result[:ind] = result_.ravel() + return result + + return cost_functional diff --git a/conmech/struct/types.py b/conmech/struct/types.py new file mode 100644 index 00000000..ac3503da --- /dev/null +++ b/conmech/struct/types.py @@ -0,0 +1,30 @@ +# CONMECH @ Jagiellonian University in Kraków +# +# Copyright (C) 2024 Piotr Bartman-Szwarc +# +# This program is free software; you can redistribute it and/or +# modify it under the terms of the GNU General Public License +# as published by the Free Software Foundation; either version 3 +# of the License, or (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, +# USA. + +from numba import types +from numba import float64 as f64 +from numba import int64 as i64 +from numba.types import Tuple + +const_i64_vec = types.Array(types.int64, 1, 'A', readonly=True) +const_i64_mat = types.Array(types.int64, 2, 'A', readonly=True) +const_f64_vec = types.Array(types.float64, 1, 'A', readonly=True) +const_f64_mat = types.Array(types.float64, 2, 'A', readonly=True) + +__all__ = ['const_i64_vec', 'const_i64_mat', 'const_f64_vec', 'const_f64_mat', 'i64', 'f64', 'Tuple'] diff --git a/examples/Bagirrov_Bartman_Ochal_2024.py b/examples/Bagirrov_Bartman_Ochal_2024.py index 2c506115..6c12bd33 100644 --- a/examples/Bagirrov_Bartman_Ochal_2024.py +++ b/examples/Bagirrov_Bartman_Ochal_2024.py @@ -16,75 +16,12 @@ from conmech.simulations.problem_solver import StaticSolver from examples.Makela_et_al_1998 import loss_value, plot_losses -mesh_density = 4 +mesh_density = 8 kN = 1000 mm = 0.001 E = 1.378e8 * kN kappa = 0.3 surface = 5 * mm * 80 * mm -k0 = 30e6 * kN * surface -k10 = 10e6 * kN * surface -k11 = 10e3 * kN * surface -k20 = 5e6 * kN * surface -k21 = 5e3 * kN * surface -k30 = 2.5e12 * kN * surface - - -class MMLV99(ContactLaw): - @staticmethod - def potential_normal_direction( - var_nu: float, static_displacement_nu: float, dt: float - ) -> float: - if var_nu <= 0: - return 0.0 - if var_nu < 0.5 * mm: - return k0 * var_nu**2 - if var_nu < 1 * mm: - return k10 * var_nu**2 + k11 * var_nu - if var_nu < 2 * mm: - return k20 * var_nu**2 + k21 * var_nu + 4 - return var_nu**4 * k30 - - @staticmethod - def potential_tangential_direction( - var_tau: float, static_displacement_tau: float, dt: float - ) -> float: - return np.log(np.sum(var_tau**2) ** 0.5 + 1) - - @staticmethod - def subderivative_normal_direction( - var_nu: float, static_displacement_nu: float, dt: float - ) -> float: - if var_nu <= 0: - return 0.0 - if var_nu < 0.5 * mm: - return k0 * var_nu * 2 - if var_nu < 1 * mm: - return k10 * (var_nu * 2) + k11 - if var_nu < 2 * mm: - return k20 * (var_nu * 2) + k21 - return var_nu**3 * 4 * k30 - - @staticmethod - def sub2derivative_normal_direction( - var_nu: float, static_displacement_nu: float, dt: float - ) -> float: - acc = 0.0 - p1 = 0.5 * mm - if var_nu <= p1: - return acc - - acc += (k0 * p1 * 2) - (k10 * (p1 * 2) + k11) - p2 = 1 * mm - if var_nu < p2: - return acc - - acc += (k10 * (p2 * 2) + k11) - (k20 * (p2 * 2) + k21) - p3 = 2 * mm - if var_nu < p3: - return acc - - return acc def make_composite_contact_law(layers, alpha, beta): @@ -160,7 +97,7 @@ class StaticSetup(StaticDisplacementProblem): elements_number: ... = (mesh_density, 8 * mesh_density) mu_coef: ... = (E * surface) / (2 * (1 + kappa)) la_coef: ... = ((E * surface) * kappa) / ((1 + kappa) * (1 - 2 * kappa)) - contact_law: ... = MMLV99 + contact_law: ... = None @staticmethod def inner_forces(x, t=None): @@ -234,7 +171,7 @@ def outer_forces(x, t=None): method=method, maxiter=100, ) - m_loss[force] = loss_value(state, validator), time.time() - start + m_loss[force] = loss_value(state, validator), runner.step_solver.computation_time path = f"{config.outputs_path}/{PREFIX}_mtd_{method}_frc_{force:.2e}" with open(path, "wb+") as output: state.body.dynamics.force.outer.source = None @@ -253,28 +190,28 @@ def outer_forces(x, t=None): # if config.show: plot_losses(path) - # for m in methods: - # for f in forces: - # path = f"{config.outputs_path}/{PREFIX}_mtd_{m}_frc_{f:.2e}" - # with open(path, "rb") as output: - # state = pickle.load(output) - # - # drawer = Drawer(state=state, config=config) - # drawer.colorful = True - # drawer.draw( - # show=config.show, - # save=config.save, - # title=f"{m}: {f}, ", - # # f"time: {runner.step_solver.last_timing}" - # ) - # x = state.body.mesh.nodes[: state.body.mesh.contact_nodes_count - 1, 0] - # u = state.displacement[: state.body.mesh.contact_nodes_count - 1, 1] - # y1 = [MMLV99().subderivative_normal_direction(-u_, 0.0, 0.0) for u_ in u] - # print(f) - # plt.plot(x, y1, label=f"{f:.2e}") - # plt.title(m) - # plt.legend() - # plt.show() + for m in methods[:]: + for f in forces[3:4]: + path = f"{config.outputs_path}/{PREFIX}_mtd_{m}_frc_{f:.2e}" + with open(path, "rb") as output: + state = pickle.load(output) + + drawer = Drawer(state=state, config=config) + drawer.colorful = True + drawer.draw( + show=config.show, + save=config.save, + title=f"{m}: {f}, ", + # f"time: {runner.step_solver.last_timing}" + ) + # x = state.body.mesh.nodes[: state.body.mesh.contact_nodes_count - 1, 0] + # u = state.displacement[: state.body.mesh.contact_nodes_count - 1, 1] + # y1 = [MMLV99().subderivative_normal_direction(-u_, 0.0, 0.0) for u_ in u] + # print(f) + # plt.plot(x, u, label=f"{f:.2e}") + # plt.title(m) + # plt.legend() + # plt.show() def composite_problem(config, layers_limits, thickness, methods, forces): @@ -290,46 +227,48 @@ def top(x): contact = make_composite_contact_law(layers_limits, alpha=bottom, beta=top) - X = np.linspace(-thickness / 4, 5 / 4 * thickness, 1000) - Y = np.empty(1000) - for i in range(1000): - Y[i] = contact.potential_normal_direction(X[i], 0.0, 0.0) - plt.plot(X, Y) - plt.show() - for i in range(1000): - Y[i] = contact.subderivative_normal_direction(X[i], 0.0, 0.0) - plt.plot(X, Y) + # X = np.linspace(-thickness / 4, 5 / 4 * thickness, 1000) + # Y = np.empty(1000) + # for i in range(1000): + # Y[i] = contact.potential_normal_direction(X[i], 0.0, 0.0) + # plt.plot(X, Y) + # plt.show() + # for i in range(1000): + # Y[i] = contact.subderivative_normal_direction(X[i], 0.0, 0.0) + # plt.plot(X, Y) + # # plt.show() + # for i in range(1000): + # Y[i] = contact.sub2derivative_normal_direction(X[i], 0.0, 0.0) + # plt.plot(X, Y) # plt.show() - for i in range(1000): - Y[i] = contact.sub2derivative_normal_direction(X[i], 0.0, 0.0) - plt.plot(X, Y) - plt.show() main(config, methods, forces, contact, prefix=str(len(layers_limits))) def survey(config): methods = ( "gradiented BFGS", - # "gradiented CG", + "gradiented CG", "BFGS", - # "CG", + "CG", "Powell", - "qsm", + # "qsm", + # "globqsm", "dc qsm", - "globqsm", "dc globqsm" )[:] forces = ( 23e3 * kN, + 24e3 * kN, 25e3 * kN, - 25.5e3 * kN, + # 25.5e3 * kN, 26e3 * kN, - 26.5e3 * kN, + # 26.5e3 * kN, 27e3 * kN, 28e3 * kN, + 29e3 * kN, 30e3 * kN, )[:] - for i in [0, 1, 2, 3, 4, 5, 6, 7, 8]: #range(1, 10 + 1, 2): + for i in [0, 1, 2, 3, 4, 5, 6, 7, 8][::2]: #range(1, 10 + 1, 2): thickness = 3 * mm layers_num = i # layers_limits = np.logspace(0, thickness, layers_num + 1) @@ -355,5 +294,5 @@ def partition(start, stop, num, p=1.0): if __name__ == "__main__": - config_ = Config(save=False, show=False, force=True).init() + config_ = Config(save=False, show=True, force=True).init() survey(config_) diff --git a/examples/Makela_et_al_1998.py b/examples/Makela_et_al_1998.py index 4f0806e3..cdeb57ad 100644 --- a/examples/Makela_et_al_1998.py +++ b/examples/Makela_et_al_1998.py @@ -316,7 +316,7 @@ def loss_value(state, runner) -> float: ( 20e3 * kN, 21e3 * kN, - 21e3 * kN, + 22e3 * kN, 23e3 * kN, 25e3 * kN, 26e3 * kN, @@ -326,5 +326,5 @@ def loss_value(state, runner) -> float: ) )[:] - methods = ("gradiented BFGS", "gradiented CG", "BFGS", "CG", "Powell", "globqsm", "qsm")[-1:] + methods = ("gradiented BFGS", "gradiented CG", "BFGS", "CG", "Powell", "globqsm", "qsm", "dc globqsm", "dc qsm")[-4:] main(Config(save=False, show=False, force=True).init(), methods, forces) From 73e0a0ed2f83a93f3af0ee20b60663330e5b98f1 Mon Sep 17 00:00:00 2001 From: Piotr Bartman-Szwarc Date: Fri, 15 Nov 2024 20:25:41 +0100 Subject: [PATCH 34/34] BBO:WIP: example 2 WIP --- examples/Bagirrov_Bartman_Ochal_2024.py | 35 +++++++++++++------------ 1 file changed, 18 insertions(+), 17 deletions(-) diff --git a/examples/Bagirrov_Bartman_Ochal_2024.py b/examples/Bagirrov_Bartman_Ochal_2024.py index 6c12bd33..0810d10b 100644 --- a/examples/Bagirrov_Bartman_Ochal_2024.py +++ b/examples/Bagirrov_Bartman_Ochal_2024.py @@ -16,12 +16,13 @@ from conmech.simulations.problem_solver import StaticSolver from examples.Makela_et_al_1998 import loss_value, plot_losses +cc = 1.5 mesh_density = 8 kN = 1000 mm = 0.001 -E = 1.378e8 * kN -kappa = 0.3 -surface = 5 * mm * 80 * mm +E = cc * 1.378e8 * kN +kappa = 0.4 +surface = 10 * mm * 100 * mm def make_composite_contact_law(layers, alpha, beta): @@ -94,7 +95,7 @@ def sub2derivative_normal_direction( @dataclass() class StaticSetup(StaticDisplacementProblem): grid_height: ... = 10 * mm - elements_number: ... = (mesh_density, 8 * mesh_density) + elements_number: ... = (mesh_density, 6 * mesh_density) mu_coef: ... = (E * surface) / (2 * (1 + kappa)) la_coef: ... = ((E * surface) * kappa) / ((1 + kappa) * (1 - 2 * kappa)) contact_law: ... = None @@ -140,7 +141,7 @@ def main(config: Config, methods, forces, contact, prefix=""): mesh_descr = RectangleMeshDescription( initial_position=None, max_element_perimeter=0.25 * 10 * mm, - scale=[8 * 10 * mm, 10 * mm], + scale=[4 * 10 * mm, 10 * mm], ) if to_simulate: @@ -223,7 +224,7 @@ def bottom(x): def top(x): # if x >= thickness: # return 0.0 - return 15.0 / mm + (x / mm) / mm + return (15.0 / mm + (x / mm) / mm) * surface / 400 / mm**2 * cc contact = make_composite_contact_law(layers_limits, alpha=bottom, beta=top) @@ -247,28 +248,28 @@ def top(x): def survey(config): methods = ( "gradiented BFGS", - "gradiented CG", + # "gradiented CG", "BFGS", - "CG", + # "CG", "Powell", - # "qsm", + "qsm", # "globqsm", - "dc qsm", + # "dc qsm", "dc globqsm" )[:] - forces = ( - 23e3 * kN, + forces = np.asarray(( + # 23e3 * kN, 24e3 * kN, - 25e3 * kN, + # 25e3 * kN, # 25.5e3 * kN, 26e3 * kN, # 26.5e3 * kN, - 27e3 * kN, + # 27e3 * kN, 28e3 * kN, - 29e3 * kN, + # 29e3 * kN, 30e3 * kN, - )[:] - for i in [0, 1, 2, 3, 4, 5, 6, 7, 8][::2]: #range(1, 10 + 1, 2): + )) * cc + for i in [0, 1, 2, 3, 4, 5, 6, 7, 8][4::2]: #range(1, 10 + 1, 2): thickness = 3 * mm layers_num = i # layers_limits = np.logspace(0, thickness, layers_num + 1)