Skip to content

Commit

Permalink
cleanup
Browse files Browse the repository at this point in the history
  • Loading branch information
Piotr Bartman-Szwarc committed Sep 24, 2024
1 parent 460b406 commit ad11179
Show file tree
Hide file tree
Showing 4 changed files with 114 additions and 109 deletions.
23 changes: 17 additions & 6 deletions conmech/plotting/drawer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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()

Expand All @@ -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
Expand Down
6 changes: 5 additions & 1 deletion conmech/solvers/optimization/optimization.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
)
Expand All @@ -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
69 changes: 21 additions & 48 deletions examples/Bagirrov_Bartman_Ochal_2024.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
"""

import pickle
import time
from dataclasses import dataclass
from matplotlib import pyplot as plt

Expand All @@ -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
Expand Down Expand Up @@ -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:
Expand All @@ -132,14 +136,17 @@ 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,
initial_displacement=setup.initial_displacement,
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
Expand All @@ -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:
Expand Down Expand Up @@ -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,
Expand All @@ -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)
125 changes: 71 additions & 54 deletions examples/Makela_et_al_1998.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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:
Expand All @@ -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
Expand All @@ -173,15 +181,22 @@ 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:
path = f"{config.outputs_path}/{PREFIX}_mtd_{m}_frc_{f:.2e}"
with open(path, "rb") as output:
state = pickle.load(output)

print("drawing")
drawer = Drawer(state=state, config=config)
drawer.colorful = True
drawer.draw(
Expand All @@ -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]
Expand All @@ -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):
Expand All @@ -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,
Expand All @@ -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)

0 comments on commit ad11179

Please sign in to comment.