From bb02f5dc820671d3cd5d37369c67df1465b63b40 Mon Sep 17 00:00:00 2001 From: Piotr Bartman Date: Sun, 18 Feb 2024 02:20:34 +1100 Subject: [PATCH] working WIP --- conmech/solvers/optimization/optimization.py | 11 ++ examples/Bagirov_Bartman_Ochal_2023.py | 119 ++++++++++++++++--- 2 files changed, 116 insertions(+), 14 deletions(-) diff --git a/conmech/solvers/optimization/optimization.py b/conmech/solvers/optimization/optimization.py index 34a5275f..a58ce532 100644 --- a/conmech/solvers/optimization/optimization.py +++ b/conmech/solvers/optimization/optimization.py @@ -77,6 +77,8 @@ def __init__( else: raise ValueError(f"Unknown statement: {statement}") + RESULTS = "" + def __str__(self): raise NotImplementedError() @@ -128,6 +130,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", @@ -150,4 +159,6 @@ def _solve_impl( solution = result.x norm = np.linalg.norm(np.subtract(solution, old_solution)) old_solution = solution.copy() + Optimization.RESULTS += \ + f'"{method}":' + str(self.loss(solution, *args)[0]) + ",\n" return solution 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)