BBO

.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

conmech/dynamics/contact/damped_normal_compliance.py

@@ -52,4 +52,13 @@ 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

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

conmech/plotting/drawer.py

@@ -5,6 +5,8 @@
 @author: Piotr Bartman
 """
 
+import time
+
 import matplotlib.pyplot as plt
 import networkx as nx
 import numpy as np
@@ -53,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()
 
@@ -76,15 +94,10 @@ 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)
-
-        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
@@ -252,7 +265,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
 

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

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):

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)

conmech/solvers/optimization/optimization.py

@@ -1,14 +1,30 @@
-"""
-Created at 18.02.2021
-"""
-
+# CONMECH @ Jagiellonian University in Kraków
+#
+# Copyright (C) 2021-2024  Piotr Bartman-Szwarc <[email protected]>
+#
+# 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
+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,
@@ -17,7 +33,29 @@
 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_subgradient, make_subgradient_dc,
+)
+
+
+QSMLM_NAMES = {
+    "quasi secant method",
+    "limited memory quasi secant method",
+    "quasi secant method limited memory",
+    "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",
+    "global quasi secant method limited memory",
+    "globqsm",
+    "globqsmlm",
+}
+GLOBAL_QSMLM_NAMES = {"dc " + name for name in GLOBAL_QSMLM_NAMES}.union(GLOBAL_QSMLM_NAMES)
 
 
 class Optimization(Solver):
@@ -44,6 +82,20 @@ def __init__(
             variable_dimension=statement.dimension_out,
             problem_dimension=statement.dimension_in,
         )
+        if hasattr(contact_law, "subderivative_normal_direction"):  # TODO
+            self.subgradient = make_subgradient(
+                normal_condition=contact_law.subderivative_normal_direction,
+            )
+        else:
+            self.subgradient = None
+        if hasattr(contact_law, "sub2derivative_normal_direction"):  # TODO
+            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
         if isinstance(statement, WaveStatement):
             if isinstance(contact_law, InteriorContactLaw):
                 self.loss = make_equation(  # TODO!
@@ -59,6 +111,7 @@ def __init__(
                 variable_dimension=statement.dimension_out,
                 problem_dimension=statement.dimension_in,
             )
+        self.minimizer = None
 
     def __str__(self):
         raise NotImplementedError()
@@ -81,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))
@@ -94,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,
@@ -104,20 +158,52 @@ 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)
+            from kosopt.qsmlm import make_minimizer
+
+            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 (
-                "quasi secant method",
-                "limited memory quasi secant method",
-                "quasi secant method limited memory",
-                "qsm",
-                "qsmlm",
+            if method.lower() in QSMLM_NAMES:
+                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, 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",
+                "dg",
             ):
                 # 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, 0, 1, 10000)
                 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
 
@@ -140,23 +226,29 @@ def constr(x):
                 solution = result.x
                 sols.append(solution)
                 loss.append(self.loss(solution, *args)[0])
-                break
             else:
+                subgrad = None
+                if method.startswith("gradiented "):
+                    subgrad = self.subgradient
+                    method = method[len("gradiented "):]
+                start = time.time()
                 result = scipy.optimize.minimize(
                     self.loss,
                     solution,
                     args=args,
                     method=method,
+                    jac=subgrad,
                     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])
-                break
 
             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]
 

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
@@ -77,8 +76,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]]

conmech/solvers/solver.py

@@ -1,7 +1,22 @@
-"""
-Created at 18.02.2021
-"""
-
+# CONMECH @ Jagiellonian University in Kraków
+#
+# Copyright (C) 2021-2024  Piotr Bartman-Szwarc <[email protected]>
+#
+# 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
@@ -48,6 +63,8 @@ def __init__(
             )
         )
 
+        self.computation_time = 0.0
+
     def __str__(self) -> str:
         raise NotImplementedError()