From ec6f3c0bdb5936bde3245b670b73ce9a7efb8d85 Mon Sep 17 00:00:00 2001 From: Piotr Date: Wed, 8 Nov 2023 10:51:21 +0100 Subject: [PATCH] BOST 2024 --- conmech/mesh/boundaries.py | 2 +- conmech/mesh/zoo/barochsoftar_2024.py | 45 +++ conmech/plotting/drawer.py | 33 +- conmech/properties/mesh_description.py | 8 + conmech/solvers/optimization/optimization.py | 47 ++- examples/BOST_2024.py | 269 ++++++++++++++ examples/Jureczka_and_Ochal_2019.py | 4 +- examples/error_estimates.py | 364 +++++-------------- 8 files changed, 489 insertions(+), 283 deletions(-) create mode 100644 conmech/mesh/zoo/barochsoftar_2024.py create mode 100644 examples/BOST_2024.py diff --git a/conmech/mesh/boundaries.py b/conmech/mesh/boundaries.py index a218a23a..4bd8041a 100644 --- a/conmech/mesh/boundaries.py +++ b/conmech/mesh/boundaries.py @@ -102,7 +102,7 @@ def contact_boundary(self): @property def contact_normals(self): - return self.surface_normals[: self.contact_nodes_count] + return self.surface_normals[: len(self.boundaries["contact"].surfaces)] @property def neumann_boundary(self): diff --git a/conmech/mesh/zoo/barochsoftar_2024.py b/conmech/mesh/zoo/barochsoftar_2024.py new file mode 100644 index 00000000..e0d849a0 --- /dev/null +++ b/conmech/mesh/zoo/barochsoftar_2024.py @@ -0,0 +1,45 @@ +# CONMECH @ Jagiellonian University in Kraków +# +# Copyright (C) 2023 Piotr Bartman +# +# 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 dmsh + +from conmech.mesh.zoo.raw_mesh import RawMesh +from conmech.properties.mesh_description import JOB2023MeshDescription + + +class BOST2023(RawMesh): + def __init__(self, mesh_descr: JOB2023MeshDescription): + # pylint: disable=no-member # for dmsh + geo = dmsh.Polygon( + [ + [0.0, 0.2], + [0.2, 0.4], + [0.6, 0.4], + [0.8, 0.2], + [0.8, 0.0], + [1.2, 0.0], + [1.2, 1.2], + [0.8, 1.2], + [0.8, 1.0], + [0.6, 0.8], + [0.2, 0.8], + [0.0, 1.0], + ] + ) + nodes, elements = dmsh.generate(geo, mesh_descr.max_element_perimeter) + super().__init__(nodes, elements) diff --git a/conmech/plotting/drawer.py b/conmech/plotting/drawer.py index 13176e2c..599a7e17 100644 --- a/conmech/plotting/drawer.py +++ b/conmech/plotting/drawer.py @@ -84,28 +84,32 @@ def draw( fig.tight_layout() plt.show() if save: - self.save_plot(save_format) + self.save_plot(save_format, name=save) def set_axes_limits(self, axes, foundation): # pylint: disable=nested-min-max if self.x_min is None: self.x_min = min( - min(self.state.body.mesh.nodes[:, 0]), min(self.state.displaced_nodes[:, 0]) + min(self.state.body.mesh.nodes[:, 0]), + min(self.state.displaced_nodes[:, 0]), ) if self.x_max is None: self.x_max = max( - max(self.state.body.mesh.nodes[:, 0]), max(self.state.displaced_nodes[:, 0]) + max(self.state.body.mesh.nodes[:, 0]), + max(self.state.displaced_nodes[:, 0]), ) dx = self.x_max - self.x_min x_margin = dx * 0.2 xlim = (self.x_min - x_margin, self.x_max + x_margin) if self.y_min is None: self.y_min = min( - min(self.state.body.mesh.nodes[:, 1]), min(self.state.displaced_nodes[:, 1]) + min(self.state.body.mesh.nodes[:, 1]), + min(self.state.displaced_nodes[:, 1]), ) if self.y_max is None: self.y_max = max( - max(self.state.body.mesh.nodes[:, 1]), max(self.state.displaced_nodes[:, 1]) + max(self.state.body.mesh.nodes[:, 1]), + max(self.state.displaced_nodes[:, 1]), ) dy = self.y_max - self.y_min y_margin = dy * 0.2 @@ -163,7 +167,10 @@ def draw_boundaries(self, axes): edges=self.mesh.contact_boundary, nodes=nodes, axes=axes, edge_color="b" ) self.draw_boundary( - edges=self.mesh.dirichlet_boundary, nodes=nodes, axes=axes, edge_color="r" + edges=self.mesh.dirichlet_boundary, + nodes=nodes, + axes=axes, + edge_color="r", ) self.draw_boundary( edges=self.mesh.neumann_boundary, nodes=nodes, axes=axes, edge_color="g" @@ -175,7 +182,7 @@ def draw_dirichlet(self, axes): dirichlet_nodes = self.state.body.mesh.dirichlet_boundary dirichlet_nodes = list(set(dirichlet_nodes.flatten())) x = self.state.displaced_nodes[dirichlet_nodes][[0, 1]] - v = np.zeros_like(x) + np.asarray([0, 1]) + v = np.zeros_like(x) + np.asarray([1, 0]) if any(v[:, 0]) or any(v[:, 1]): # to avoid warning axes.quiver( x[:, 0], @@ -184,7 +191,7 @@ def draw_dirichlet(self, axes): v[:, 1], angles="xy", scale_units="xy", - scale=2.5, + scale=8.0, headlength=10, headaxislength=10, headwidth=10, @@ -245,7 +252,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 name else config.timestamp + name = name if isinstance(name, str) else config.timestamp path = f"{directory}/{name}.{format_}" return path @@ -261,7 +268,9 @@ def save_plot(self, format_, name=None): ) plt.close() - def draw_boundary(self, edges, nodes, axes, label="", node_color="k", edge_color="k"): + def draw_boundary( + self, edges, nodes, axes, label="", node_color="k", edge_color="k" + ): graph = nx.Graph() for edge in edges: graph.add_edge(edge[0], edge[1]) @@ -303,6 +312,8 @@ def draw_field(self, field, v_min, v_max, axes, fig): # from mpl_toolkits.axes_grid1 import make_axes_locatable # divider = make_axes_locatable(axes) # cax = divider.append_axes("bottom", size="5%", pad=0.15) - sm = plt.cm.ScalarMappable(cmap=self.cmap, norm=plt.Normalize(vmin=v_min, vmax=v_max)) + sm = plt.cm.ScalarMappable( + cmap=self.cmap, norm=plt.Normalize(vmin=v_min, vmax=v_max) + ) sm.set_array([]) fig.colorbar(sm, orientation="horizontal", label=self.field_label, ax=axes) diff --git a/conmech/properties/mesh_description.py b/conmech/properties/mesh_description.py index c8afb649..f9877820 100644 --- a/conmech/properties/mesh_description.py +++ b/conmech/properties/mesh_description.py @@ -89,6 +89,14 @@ def build(self): return Barboteu2008(self) +@dataclass +class BOST2023MeshDescription(GeneratedMeshDescription): + def build(self): + from conmech.mesh.zoo.barochsoftar_2024 import BOST2023 + + return BOST2023(self) + + @dataclass class JOB2023MeshDescription(GeneratedMeshDescription): def build(self): diff --git a/conmech/solvers/optimization/optimization.py b/conmech/solvers/optimization/optimization.py index 4ad30c12..4c80716b 100644 --- a/conmech/solvers/optimization/optimization.py +++ b/conmech/solvers/optimization/optimization.py @@ -93,8 +93,8 @@ def _solve_impl( displacement = np.squeeze(displacement.copy().reshape(1, -1)) old_solution = np.squeeze(initial_guess.copy().reshape(1, -1)) disp = kwargs.get("disp", False) - maxiter = kwargs.get("maxiter", int(len(initial_guess) * 1e9)) - tol = kwargs.get("tol", 1e-12) + maxiter = kwargs.get("maxiter", int(len(initial_guess) * 1e12)) + tol = kwargs.get("tol", 1e-18) args = ( self.body.mesh.nodes, self.body.mesh.contact_boundary, @@ -105,6 +105,13 @@ def _solve_impl( self.time_step, ) + loss = [] + sols = [] + sols.append(solution) + loss.append(self.loss(solution, *args)[0]) + print("pre", method, self.loss(solution, *args)) + + while norm >= fixed_point_abs_tol: if method.lower() in ( # TODO "quasi secant method", @@ -117,6 +124,34 @@ def _solve_impl( from kosopt import qsmlm solution = qsmlm.minimize(self.loss, solution, args=args, maxiter=maxiter) + sols[self.loss(solution, *args)] = solution + elif method.lower() == 'constrained': + import numba + contact_nodes_count = self.body.mesh.boundaries.contact_nodes_count + GAP = 0.0 + @numba.njit() + def constr(x): + offset = len(x) // 2 + t = x[offset:offset+contact_nodes_count] + return np.min(t + GAP) + + maxiter = kwargs.get("maxiter", int(len(initial_guess) * 1e9)) + tol = kwargs.get("tol", 1e-12) + result = scipy.optimize.minimize( + self.loss, + solution, + args=args, + options={"disp": disp, "maxiter": maxiter}, + tol=tol, + constraints=({'type': 'ineq', + 'fun': constr + }) + ) + solution = result.x + sols.append(solution) + loss.append(self.loss(solution, *args)[0]) + print(method, self.loss(solution, *args)) + break else: result = scipy.optimize.minimize( self.loss, @@ -127,6 +162,14 @@ def _solve_impl( tol=tol, ) solution = result.x + sols.append(solution.copy()) + loss.append(self.loss(solution, *args)[0]) + print(method, self.loss(solution, *args)) + break + norm = np.linalg.norm(np.subtract(solution, old_solution)) old_solution = solution.copy() + min_index = loss.index(np.min(loss)) + solution = sols[min_index] + print(loss, "final:", self.loss(solution, *args)) return solution diff --git a/examples/BOST_2024.py b/examples/BOST_2024.py new file mode 100644 index 00000000..068b1586 --- /dev/null +++ b/examples/BOST_2024.py @@ -0,0 +1,269 @@ +# CONMECH @ Jagiellonian University in Kraków +# +# Copyright (C) 2023 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 +from typing import Iterable, List, Optional + +import matplotlib.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.scenarios.problems import ContactLaw, StaticDisplacementProblem +from conmech.simulations.problem_solver import StaticSolver as StaticProblemSolver +from conmech.properties.mesh_description import BOST2023MeshDescription +from examples.error_estimates import error_estimates + +GAP = 0.0 +E = 12000 +kappa = 0.42 + + +class BOST23(ContactLaw): + @staticmethod + def potential_normal_direction(u_nu: float) -> float: + raise NotImplementedError() + + @staticmethod + def potential_tangential_direction(u_tau: np.ndarray) -> float: + return np.sum(u_tau * u_tau) ** 0.5 + + @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): + mu_coef: ... = E / (1 + kappa) + la_coef: ... = E * kappa / ((1 + kappa) * (1 - 2 * kappa)) + contact_law: ... = BOST23 + + @staticmethod + def inner_forces(x, t=None): + return np.array([-0.2, -0.8]) * 1e3 + + @staticmethod + def outer_forces(x, t=None): + return np.array([0, 0]) + + @staticmethod + def friction_bound(u_nu: float) -> float: + if u_nu < 0: + return 0 + return u_nu + + boundaries: ... = BoundariesDescription( + contact=lambda x: x[1] == 0, dirichlet=lambda x: x[0] == 0 + ) + + +def prepare_setup(ig, setup): + if ig == "inf": + def potential_normal_direction(u_nu: float) -> float: + return 0 + + kwargs = {"method": "constrained"} + else: + def potential_normal_direction(u_nu: float) -> float: + u_nu -= GAP + # EXAMPLE 10 + a = 0.1 + b = 0.1 + if u_nu <= 0: + result = 0.0 + elif u_nu < b: + result = (a + np.exp(-b)) / (2 * b) * u_nu ** 2 + else: + result = a * u_nu - np.exp(- u_nu) + ((b + 2) * np.exp(-b) - a * b) / 2 + return ig * result + + kwargs = {"method": "POWELL"} + + setup.contact_law.potential_normal_direction = potential_normal_direction + return kwargs + + +def main(config: Config, igs: List[int], highlighted: Iterable): + """ + Entrypoint to example. + + To see result of simulation you need to call from python `main(Config().init())`. + """ + PREFIX = "BOST_POWELL" + to_simulate = [] + for ig in igs: + try: + with open(f"{config.outputs_path}/{PREFIX}_ig_{ig}", "rb") as output: + _ = pickle.load(output) + except IOError as ioe: + print(ioe) + to_simulate.append(ig) + if config.force: + to_simulate = igs + + mesh_descr = BOST2023MeshDescription( + initial_position=None, + max_element_perimeter=1 / 32 + ) + + if to_simulate: + print("Simulating...") + print(to_simulate) + setup = StaticSetup(mesh_descr) + + if to_simulate[0] == igs[0]: + initial_displacement = setup.initial_displacement + else: + ig_prev = igs[igs.index(to_simulate[0]) - 1] + with open(f"{config.outputs_path}/{PREFIX}_ig_{ig_prev}", "rb") as output: + state = pickle.load(output) + initial_displacement = lambda _: state.displacement.copy() + + for i, ig in enumerate(to_simulate): + kwargs = prepare_setup(ig, setup) + if config.test: + setup.elements_number = (2, 4) + runner = StaticProblemSolver(setup, "schur") + + state = runner.solve( + verbose=True, + fixed_point_abs_tol=0.001, + initial_displacement=initial_displacement, + **kwargs, + ) + + if i + 1 < len(to_simulate): + ig_prev = igs[igs.index(to_simulate[i + 1]) - 1] + if ig == ig_prev: + initial_displacement = setup.initial_displacement + else: + with open(f"{config.outputs_path}/{PREFIX}_ig_{ig_prev}", "rb") as output: + state = pickle.load(output) + initial_displacement = lambda _: state.displacement.copy() + + state.displaced_nodes[: state.body.mesh.nodes_count, :] = ( + state.body.mesh.nodes[: state.body.mesh.nodes_count, :] + state.displacement[:, :] + ) + with open(f"{config.outputs_path}/{PREFIX}_ig_{ig}", "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(f"Plotting {igs=}") + for i, ig in enumerate(igs): + with open(f"{config.outputs_path}/{PREFIX}_ig_{ig}", "rb") as output: + state = pickle.load(output) + + state.displaced_nodes[:, 1] += GAP + state.body.mesh.nodes[:, 1] += GAP + + if ig in highlighted: + drawer = Drawer(state=state, config=config) + drawer.colorful = False + drawer.outer_forces_scale = 1 + if ig == 0: + title = r"$\lambda^{-1} = 0.000$" + else: + value = f"{float(np.log10(ig)):.3f}" if isinstance(ig, int) else str(ig) + title = r"$\log_{10}\lambda^{-1} = $" + value + if config.save: + save = str(ig) + else: + save = False + drawer.draw(show=config.show, save=save, title=title) + + print("Error estimate") + errors = error_estimates( + f"{config.outputs_path}/{PREFIX}_ig_inf", + *[f"{config.outputs_path}/{PREFIX}_ig_{ig}" for ig in igs]) + X = [ig for ig in igs if not isinstance(ig, str)] + Y = list(errors.values())[:-1] + Y = [v[0] for v in Y] + Y = - np.asarray(Y) + plot_errors(X, Y, highlighted_id=None, save=f"{config.outputs_path}/convergence.pdf") + + +def plot_errors(X, Y, highlighted_id, save: Optional[str] = None): + plt.plot(X[:], Y[:], marker='o', color='gray') + if highlighted_id is not None: + plt.plot(X[highlighted_id], Y[highlighted_id], 'ro', color='black') + plt.loglog() + plt.grid(True, which='major', linestyle='--', linewidth=0.5) + plt.xlabel(r"$\lambda^{-1}_{\,\, n}$") + plt.ylabel(r"$||\mathbf{u}^h_n - \mathbf{u}||$") + if save is None: + plt.show() + else: + plt.savefig(save, format="pdf") + + +if __name__ == "__main__": + + # results from the paper + highlighted = (1, 1024, 1216, 2048, 5931641, "inf") + highlighted_id = np.asarray((1, 20, 24, 26, -3)) + + X = np.asarray([0, 1, 2, 4, 5, 8, 11, 16, 22, 32, 45, 64, 90, 128, 181, 256, 362, 512, 724, 896, 1024, + 1088, 1120, 1152, 1216, 1448, 2048, 2896, 4096, 5792, 8192, 11585, 16384, 23170, 32768, + 46340, 65536, 92681, 131072, 185363, 262144, 370727, 524288, 741455, 1048576, 1482910, + 2097152, 2965820, 4194304, 5931641, 8388608, 11863283]) + Y = np.asarray([1.957688239442841, 1.955823800273701, 1.9550737995244916, 1.9535713010795799, + 1.9528202182876828, 1.9505617554633947, 1.9483005592237832, 1.9445216877172742, + 1.939976356940407, 1.9323603026049123, 1.9224029295661968, 1.9077076426739854, + 1.887348094110511, 1.8570040907627985, 1.8135169713254757, 1.7493452687885334, + 1.6528260472436915, 1.5016916054916745, 1.247485589366432, 0.9833463022514393, + 0.7219304912004216, 0.5510492874873285, 0.44668589874753256, + 0.32141797154720314, 0.2235180524765496, 0.18110766227417885, + 0.13217099205041302, 0.09590338285207281, 0.06933759603103229, + 0.050067258317548866, 0.036147047768223, 0.02613519925974203, + 0.018945653731469527, 0.013792699500121302, 0.01010466081623074, + 0.0074703782466904725, 0.0055919695751946025, 0.004256100793696483, + 0.003309243652330325, 0.002641490469903024, 0.002173650857273765, + 0.0018483357409417435, 0.0016235440683907585, 0.0014687173119601173, + 0.0013620445313966081, 0.0012883643852057848, 0.0012373372831507256, + 0.0012019371089261616, 0.0011773413958895045, 0.0011602208824840082, + 0.0011019352332271976, 0.001139913361649654] + ) + + plot_errors(X, Y, highlighted_id, save='convergence.pdf') + + show = True + + igs = {int(2**(i/2)) for i in range(0, 48)} + eigs = {1024-128, 1024+64, 1024+96, 1024+128, 1024+128+64,} + igs.update(eigs) + igs = [0] + sorted(list(igs)) + ["inf"]# + igs = igs[:] + main( + Config(save=not show, show=show, force=False).init(), + igs, highlighted + ) diff --git a/examples/Jureczka_and_Ochal_2019.py b/examples/Jureczka_and_Ochal_2019.py index 800e79cf..80dc7760 100644 --- a/examples/Jureczka_and_Ochal_2019.py +++ b/examples/Jureczka_and_Ochal_2019.py @@ -50,11 +50,11 @@ class StaticSetup(StaticDisplacementProblem): @staticmethod def inner_forces(x, t=None): - return np.array([-1.2, -0.8]) + return np.array([-1.2, -0.9]) @staticmethod def outer_forces(x, t=None): - return np.array([0, 0]) + return np.array([0.5, 0.5]) @staticmethod def friction_bound(u_nu: float) -> float: diff --git a/examples/error_estimates.py b/examples/error_estimates.py index 774d0466..bfe462f5 100644 --- a/examples/error_estimates.py +++ b/examples/error_estimates.py @@ -21,7 +21,7 @@ def compare(ref: TemperatureState, sol: TemperatureState): soltri = tri.Triangulation(x, y, triangles=sol.body.mesh.elements) u1hi = tri.LinearTriInterpolator(soltri, sol.displacement[:, 0]) u2hi = tri.LinearTriInterpolator(soltri, sol.displacement[:, 1]) - thi = tri.LinearTriInterpolator(soltri, sol.temperature) + # thi = tri.LinearTriInterpolator(soltri, sol.temperature) for element in ref.body.mesh.elements: x0 = ref.body.mesh.nodes[element[0]] @@ -35,17 +35,19 @@ def compare(ref: TemperatureState, sol: TemperatureState): u2_1 = ref.displacement[element[1], 1] u2_2 = ref.displacement[element[2], 1] u2 = u2_0 + u2_1 + u2_2 - t0 = ref.temperature[element[0]] - t1 = ref.temperature[element[1]] - t2 = ref.temperature[element[2]] - t = t0 + t1 + t2 + # t0 = ref.temperature[element[0]] + # t1 = ref.temperature[element[1]] + # t2 = ref.temperature[element[2]] + # t = t0 + t1 + t2 u1dx, u1dy = calculate_dx_dy(x0, u1_0, x1, u1_1, x2, u1_2) u2dx, u2dy = calculate_dx_dy(x0, u2_0, x1, u2_1, x2, u2_2) - tdx, tdy = calculate_dx_dy(x0, t0, x1, t1, x2, t2) - u1hdx, u1hdy = calculate_dx_dy(x0, u1hi(*x0), x1, u1hi(*x1), x2, u1hi(*x2)) - u2hdx, u2hdy = calculate_dx_dy(x0, u2hi(*x0), x1, u2hi(*x1), x2, u2hi(*x2)) - thdx, thdy = calculate_dx_dy(x0, thi(*x0), x1, thi(*x1), x2, thi(*x2)) + # tdx, tdy = calculate_dx_dy(x0, t0, x1, t1, x2, t2) + u1hdx, u1hdy = calculate_dx_dy( + x0, u1hi(*x0).compressed(), x1, u1hi(*x1).compressed(), x2, u1hi(*x2).compressed()) + u2hdx, u2hdy = calculate_dx_dy( + x0, u2hi(*x0).compressed(), x1, u2hi(*x1).compressed(), x2, u2hi(*x2).compressed()) + # thdx, thdy = calculate_dx_dy(x0, thi(*x0), x1, thi(*x1), x2, thi(*x2)) u1h = u1hi(*x0) + u1hi(*x1) + u1hi(*x2) u2h = u2hi(*x0) + u2hi(*x1) + u2hi(*x2) @@ -57,8 +59,8 @@ def compare(ref: TemperatureState, sol: TemperatureState): + (u1dy - u1hdy + u2dx - u2hdx) ** 2 ) ** 0.5 - th = thi(*x0) + thi(*x1) + thi(*x2) - tt += ((t - th) ** 2 + (tdx - thdx) ** 2 + (tdy - thdy) ** 2) ** 0.5 + # th = thi(*x0) + thi(*x1) + thi(*x2) + # tt += ((t - th) ** 2 + (tdx - thdx) ** 2 + (tdy - thdy) ** 2) ** 0.5 return ut, tt @@ -78,268 +80,96 @@ def calculate_dx_dy(x0, u0, x1, u1, x2, u2): return dx, dy -if __name__ == "__main__": +def error_estimates(ref, *args): reference_k_h = (9, 6) - denominator = 2 ** reference_k_h[0] * 2 ** reference_k_h[1] * 4 + denominator = 1 # TODO 2 ** reference_k_h[0] * 2 ** reference_k_h[1] * 4 - T = 1 - kn = 10 - hn = 6 - ks = range(kn) - hs = range(hn) + # T = 1 + # kn = 10 + # hn = 6 + # ks = range(kn) + # hs = range(hn) - # with open(f'output/temp/k_{reference_k_h[0]}_h_{reference_k_h[1]}', 'rb') as output: - # reference = pickle.load(output) - # - # ue = np.empty((kn, hn)) + with open(ref, 'rb') as output: + reference = pickle.load(output) + denominator = len(reference.body.mesh.elements) + print(f"{denominator=}") + + ue = {} # te = np.empty((kn, hn)) - # for h in hs: - # for k in ks: - # with open(f'output/temp/k_{k}_h_{h}', 'rb') as output: - # solution = pickle.load(output) - # u, t = compare(reference, solution) - # ue[k, h] = u / denominator - # te[k, h] = t / denominator - # print(k, h, u, t) + for arg in args: + with open(arg, 'rb') as output: + solution = pickle.load(output) + u, t = compare(reference, solution) + ue[arg] = u / denominator + # te[k, h] = t / denominator + # print(k, h, u, t) + return ue # print(repr(ue)) - # print(repr(te)) - # ue = np.asarray([[3.51064292e-03, 6.60816772e-04, 1.57858401e-04, 8.49138212e-05, - # 7.25934268e-05, 7.06925935e-05], - # [3.61655733e-03, 6.48988624e-04, 1.15492804e-04, 3.93501894e-05, - # 2.77189199e-05, 2.61284821e-05], - # [3.76727797e-03, 6.54611297e-04, 1.12705825e-04, 2.88279473e-05, - # 1.84614324e-05, 1.72210597e-05], - # [3.67165437e-03, 6.50281211e-04, 1.04942153e-04, 2.40933567e-05, - # 1.22119143e-05, 9.96705359e-06], - # [3.55872410e-03, 6.90665486e-04, 9.86944361e-05, 1.94383337e-05, - # 6.82798321e-06, 4.58554490e-06], - # [3.38251554e-03, 6.68536915e-04, 9.80099335e-05, 1.72279032e-05, - # 5.17652041e-06, 1.90904509e-06], - # [3.63972474e-03, 6.55831166e-04, 9.67532875e-05, 1.64753487e-05, - # 3.39027433e-06, 1.14619503e-06], - # [3.84269955e-03, 6.57062689e-04, 9.43591112e-05, 1.60484103e-05, - # 3.06911182e-06, 6.50277890e-07], - # [3.89351253e-03, 6.53322099e-04, 9.78647821e-05, 2.84040210e-05, - # 3.21183711e-06, 6.66848485e-07], - # [3.56222413e-03, 7.58684791e-04, 1.01876870e-04, 1.69102222e-05, - # 3.16765328e-06, 6.76218102e-07]]) - # te = np.asarray([[2.51983950e-04, 3.65203893e-05, 1.40350369e-05, 9.03639485e-06, - # 7.95321197e-06, 7.68786449e-06], - # [2.58693772e-04, 2.74905815e-05, 7.45569008e-06, 3.86424203e-06, - # 3.13708124e-06, 2.96597709e-06], - # [2.71314633e-04, 2.36231095e-05, 4.67394194e-06, 1.76615566e-06, - # 1.24542801e-06, 1.13771908e-06], - # [2.62131942e-04, 2.20884726e-05, 3.55819379e-06, 9.05536091e-07, - # 4.66154672e-07, 3.82321131e-07], - # [2.53106093e-04, 2.17287469e-05, 3.17266640e-06, 6.22052400e-07, - # 2.10371400e-07, 1.25786283e-07], - # [2.37256953e-04, 2.13680291e-05, 3.10738659e-06, 5.26606195e-07, - # 1.20051134e-07, 5.10547779e-08], - # [2.56370262e-04, 2.11524557e-05, 3.00687181e-06, 4.92886048e-07, - # 1.06324965e-07, 2.91076083e-08], - # [2.79046413e-04, 2.09837016e-05, 3.00190538e-06, 4.83016121e-07, - # 9.65938350e-08, 2.08012631e-08], - # [2.83269851e-04, 2.13542387e-05, 3.13429246e-06, 6.21141133e-07, - # 9.40967476e-08, 1.82546450e-08], - # [2.54791829e-04, 2.08020592e-05, 2.96826861e-06, 4.77383581e-07, - # 9.37414078e-08, 1.81747962e-08]]) - # array([[3.14681433e-03, 5.78708271e-04, 1.42346139e-04, 7.61790675e-05, - # 6.42927759e-05, 6.19299357e-05], - # [3.25632082e-03, 5.39096000e-04, 9.65981575e-05, 3.50823920e-05, - # 2.57918955e-05, 2.45186524e-05], - # [3.41123943e-03, 5.28468199e-04, 9.15348410e-05, 2.56499495e-05, - # 1.81253651e-05, 1.72970692e-05], - # [3.31525800e-03, 5.17078003e-04, 8.21106003e-05, 2.06276332e-05, - # 1.18540371e-05, 1.01068609e-05], - # [3.19620516e-03, 5.15337705e-04, 7.44393552e-05, 1.56228857e-05, - # 6.31159301e-06, 4.60707382e-06], - # [3.01144473e-03, 5.09017138e-04, 7.29054520e-05, 1.30637947e-05, - # 4.26555682e-06, 1.84300202e-06], - # [3.26973059e-03, 5.08055560e-04, 7.11474149e-05, 1.21325368e-05, - # 2.61712020e-06, 1.03775258e-06], - # [3.48839816e-03, 5.15453847e-04, 6.74653739e-05, 1.14800101e-05, - # 2.28290731e-06, 5.13143518e-07], - # [3.53998890e-03, 5.05289847e-04, 7.01592614e-05, 1.97090876e-05, - # 2.37363292e-06, 5.21670964e-07], - # [3.20002421e-03, 5.59609645e-04, 7.12836450e-05, 1.17558040e-05, - # 2.35212189e-06, 5.25712091e-07]]) - # array([[2.51983950e-04, 3.65203893e-05, 1.40350369e-05, 9.03639485e-06, - # 7.95321197e-06, 7.68786449e-06], - # [2.58693772e-04, 2.74905815e-05, 7.45569008e-06, 3.86424203e-06, - # 3.13708124e-06, 2.96597709e-06], - # [2.71314633e-04, 2.36231095e-05, 4.67394194e-06, 1.76615566e-06, - # 1.24542801e-06, 1.13771908e-06], - # [2.62131942e-04, 2.20884726e-05, 3.55819379e-06, 9.05536091e-07, - # 4.66154672e-07, 3.82321131e-07], - # [2.53106093e-04, 2.17287469e-05, 3.17266640e-06, 6.22052400e-07, - # 2.10371400e-07, 1.25786283e-07], - # [2.37256953e-04, 2.13680291e-05, 3.10738659e-06, 5.26606195e-07, - # 1.20051134e-07, 5.10547779e-08], - # [2.56370262e-04, 2.11524557e-05, 3.00687181e-06, 4.92886048e-07, - # 1.06324965e-07, 2.91076083e-08], - # [2.79046413e-04, 2.09837016e-05, 3.00190538e-06, 4.83016121e-07, - # 9.65938350e-08, 2.08012631e-08], - # [2.83269851e-04, 2.13542387e-05, 3.13429246e-06, 6.21141133e-07, - # 9.40967476e-08, 1.82546450e-08], - # [2.54791829e-04, 2.08020592e-05, 2.96826861e-06, 4.77383581e-07, - # 9.37414078e-08, 1.81747962e-08]]) - ue = np.asarray( - [ - [0.02243367, 0.00945494, 0.0047539, 0.00355776, 0.00331022, 0.00324383], - [0.02291673, 0.00892371, 0.00371959, 0.00238767, 0.00211572, 0.00208179], - [0.02358811, 0.008724, 0.0035803, 0.00198949, 0.00178018, 0.00175999], - [0.02318723, 0.0086035, 0.00331536, 0.00177648, 0.00144094, 0.0013479], - [0.02264709, 0.00844671, 0.00308759, 0.0014665, 0.0010232, 0.00091025], - [0.02179787, 0.00843854, 0.00304546, 0.00128954, 0.0008246, 0.00056634], - [0.02290356, 0.00848389, 0.00299967, 0.00122134, 0.00058402, 0.00041183], - [0.023919, 0.00860834, 0.00283703, 0.00116145, 0.00051856, 0.00026714], - [0.02413345, 0.00849438, 0.00292525, 0.00167922, 0.00054691, 0.00027684], - [0.0226687, 0.00899269, 0.00298859, 0.00119544, 0.00054671, 0.00027651], - ] - ) - te = np.asarray( - [ - [ - 6.82540994e-03, - 2.56319187e-03, - 1.55330144e-03, - 1.21156666e-03, - 1.12150221e-03, - 1.09839411e-03, - ], - [ - 6.91144671e-03, - 2.20870880e-03, - 1.11799969e-03, - 7.59382694e-04, - 6.59473669e-04, - 6.33192738e-04, - ], - [ - 7.08020402e-03, - 2.04717024e-03, - 8.88736589e-04, - 5.11961797e-04, - 4.05334124e-04, - 3.77880307e-04, - ], - [ - 6.94868359e-03, - 1.98380168e-03, - 7.76518636e-04, - 3.73911777e-04, - 2.53195833e-04, - 2.20284279e-04, - ], - [ - 6.81657653e-03, - 1.97525209e-03, - 7.31400392e-04, - 3.12309945e-04, - 1.74347377e-04, - 1.29991907e-04, - ], - [ - 6.58256367e-03, - 1.95825699e-03, - 7.23236071e-04, - 2.85801473e-04, - 1.31119651e-04, - 8.54069481e-05, - ], - [ - 6.86190756e-03, - 1.94632624e-03, - 7.10184903e-04, - 2.74813899e-04, - 1.23550007e-04, - 6.53353204e-05, - ], - [ - 7.18011214e-03, - 1.93730438e-03, - 7.08811210e-04, - 2.70437905e-04, - 1.16131191e-04, - 5.46813675e-05, - ], - [ - 7.23629516e-03, - 1.95936459e-03, - 7.25445385e-04, - 3.02455407e-04, - 1.13763883e-04, - 5.06783627e-05, - ], - [ - 6.84668307e-03, - 1.92875550e-03, - 7.02538581e-04, - 2.66845080e-04, - 1.13258321e-04, - 5.03505677e-05, - ], - ] - ) + # print(repr(te)) - h_ticks = [1 / 2**h for h in hs] - plt.style.use("seaborn") - sns.set(rc={"axes.facecolor": "#E6EDF4"}) + # h_ticks = [1 / 2 ** h for h in hs] + # plt.style.use("seaborn") + # sns.set(rc={"axes.facecolor": "#E6EDF4"}) + # + # plt.xscale("log") + # plt.yscale("log") + # plt.title("Velocity error") + # plt.xlabel("spatial step: h") + # plt.ylabel("error: $||u-u^h||_V$") + # optimal = 2 * ue[ks[-1], 1] * np.asarray(h_ticks) ** 1 + # plt.plot(h_ticks, optimal, color="silver", linewidth=4.0, label=None) + # colors = pl.cm.jet(np.linspace(0, 1, len(ks))) + # for i, k in enumerate(ks): + # plt.plot( + # h_ticks, ue[k, hs[0]: hs[-1] + 1], "s-", label="$2^{-" + str(k) + "}$", color=colors[i] + # ) + # plt.xticks(h_ticks, h_ticks) + # plt.legend(title="time step: k") + # margin = 0.15 + # plt.xlim(max(h_ticks) * (1 + margin), min(h_ticks) * (1 - margin)) + # plt.savefig( + # "output/error/displacement.pdf", + # transparent=False, + # bbox_inches="tight", + # format="pdf", + # pad_inches=0.1, + # dpi=800, + # ) + # plt.show() + + # plt.xscale("log") + # plt.yscale("log") + # plt.title("Temperature error") + # plt.xlabel("spatial step: h") + # plt.ylabel(r"error: $||\theta-\theta^h||_E$") + # plt.plot( + # h_ticks, + # (2 * te[ks[-1], 1] * np.asarray(h_ticks) ** 1), + # color="silver", + # linewidth=4.0, + # label=None, + # ) + # colors = pl.cm.jet(np.linspace(0, 1, len(ks))) + # for i, k in enumerate(ks): + # plt.plot( + # h_ticks, te[k, hs[0]: hs[-1] + 1], "s-", label="$2^{-" + str(k) + "}$", color=colors[i] + # ) + # plt.xticks(h_ticks, h_ticks) + # plt.legend(title="time step: k") + # margin = 0.15 + # plt.xlim(max(h_ticks) * (1 + margin), min(h_ticks) * (1 - margin)) + # plt.savefig( + # "output/error/temperature.pdf", + # transparent=False, + # bbox_inches="tight", + # format="pdf", + # pad_inches=0.1, + # dpi=800, + # ) + # plt.show() - plt.xscale("log") - plt.yscale("log") - plt.title("Velocity error") - plt.xlabel("spatial step: h") - plt.ylabel("error: $||v-v^h||_V$") - optimal = 2 * ue[ks[-1], 1] * np.asarray(h_ticks) ** 1 - plt.plot(h_ticks, optimal, color="silver", linewidth=4.0, label=None) - colors = pl.cm.jet(np.linspace(0, 1, len(ks))) - for i, k in enumerate(ks): - plt.plot( - h_ticks, ue[k, hs[0] : hs[-1] + 1], "s-", label="$2^{-" + str(k) + "}$", color=colors[i] - ) - plt.xticks(h_ticks, h_ticks) - plt.legend(title="time step: k") - margin = 0.15 - plt.xlim(max(h_ticks) * (1 + margin), min(h_ticks) * (1 - margin)) - plt.savefig( - "output/error/velocity.pdf", - transparent=False, - bbox_inches="tight", - format="pdf", - pad_inches=0.1, - dpi=800, - ) - plt.show() - plt.xscale("log") - plt.yscale("log") - plt.title("Temperature error") - plt.xlabel("spatial step: h") - plt.ylabel(r"error: $||\theta-\theta^h||_E$") - plt.plot( - h_ticks, - (2 * te[ks[-1], 1] * np.asarray(h_ticks) ** 1), - color="silver", - linewidth=4.0, - label=None, - ) - colors = pl.cm.jet(np.linspace(0, 1, len(ks))) - for i, k in enumerate(ks): - plt.plot( - h_ticks, te[k, hs[0] : hs[-1] + 1], "s-", label="$2^{-" + str(k) + "}$", color=colors[i] - ) - plt.xticks(h_ticks, h_ticks) - plt.legend(title="time step: k") - margin = 0.15 - plt.xlim(max(h_ticks) * (1 + margin), min(h_ticks) * (1 - margin)) - plt.savefig( - "output/error/temperature.pdf", - transparent=False, - bbox_inches="tight", - format="pdf", - pad_inches=0.1, - dpi=800, - ) - plt.show() +if __name__ == "__main__": + pass