From ed4d7a9454beef54d5b2396959c1268ac2ba9a95 Mon Sep 17 00:00:00 2001 From: Pierric Mora Date: Wed, 3 Apr 2024 00:23:41 +0200 Subject: [PATCH 01/12] add light build option --- docsrc/conf.py | 4 ++++ setup.py | 7 ++++--- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/docsrc/conf.py b/docsrc/conf.py index 70f9bcd..be89f1e 100644 --- a/docsrc/conf.py +++ b/docsrc/conf.py @@ -51,6 +51,10 @@ nbsphinx_custom_formats = { '.py': ['jupytext.reads', {'fmt': 'py:light'}], } +if os.environ.get('ELASTODYNAMICSX_DOCS_LIGHT', 'FALSE').upper() == 'TRUE': + nbsphinx_execute = 'never' +else: + nbsphinx_execute = 'auto' copybutton_exclude = '.linenos, .gp' diff --git a/setup.py b/setup.py index cfee93f..bae009c 100644 --- a/setup.py +++ b/setup.py @@ -1,6 +1,6 @@ import setuptools import subprocess -# import os +import os with open("Readme.md", "r") as fh: long_description = fh.read() @@ -8,16 +8,17 @@ class MakeTheDoc(setuptools.Command): description = "Generate Documentation Pages using Sphinx" - user_options = [] + user_options = [('light', 'l', 'light build')] def initialize_options(self): - pass + self.light = None def finalize_options(self): pass def run(self): """The command to run when users invoke python setup.py doc""" + os.environ['ELASTODYNAMICSX_DOCS_LIGHT'] = 'TRUE' if self.light else 'FALSE' subprocess.run( ['sphinx-build docsrc docs'], shell=True) subprocess.run( From 112d7ad01f5bcd81c2a4655135ad417b8d5f3e9b Mon Sep 17 00:00:00 2001 From: Pierric Mora Date: Fri, 5 Apr 2024 00:44:59 +0200 Subject: [PATCH 02/12] select thumbnail img --- demo/freq_2D-SH_FullSpace.py | 1 + 1 file changed, 1 insertion(+) diff --git a/demo/freq_2D-SH_FullSpace.py b/demo/freq_2D-SH_FullSpace.py index 433075f..24b2936 100644 --- a/demo/freq_2D-SH_FullSpace.py +++ b/demo/freq_2D-SH_FullSpace.py @@ -1,6 +1,7 @@ # --- # jupyter: # jupytext: +# cell_metadata_json: true # text_representation: # extension: .py # format_name: light From b91c3cef3bb6234d579534c650bdcaae68621532 Mon Sep 17 00:00:00 2001 From: Pierric Mora Date: Fri, 5 Apr 2024 00:45:43 +0200 Subject: [PATCH 03/12] (some) jupyter execute instead of text --- docsrc/gettingstarted/overview.rst | 58 ++++++++++++++++++++++++------ 1 file changed, 47 insertions(+), 11 deletions(-) diff --git a/docsrc/gettingstarted/overview.rst b/docsrc/gettingstarted/overview.rst index fbf9dc0..463011a 100644 --- a/docsrc/gettingstarted/overview.rst +++ b/docsrc/gettingstarted/overview.rst @@ -18,9 +18,40 @@ Build problems -------------- Using the ``elastodynamicsx.pde`` package: +.. jupyter-execute:: + :hide-code: + + import numpy as np + + from dolfinx import mesh, fem, default_scalar_type + from mpi4py import MPI + + degElement = 1 + length, height = 1, 1 + Nx, Ny = 10//degElement, 10//degElement + + # create the mesh + extent = [[0., 0.], [length, height]] + domain = mesh.create_rectangle(MPI.COMM_WORLD, extent, [Nx, Ny], mesh.CellType.triangle) + + # create the function space + V = fem.FunctionSpace(domain, ("Lagrange", degElement, (domain.geometry.dim,))) + + from elastodynamicsx.utils import make_facet_tags, make_cell_tags + # define some tags + tag_top = 1 + boundaries = [(tag_top, lambda x: np.isclose(x[1], height)),] + subdomains = [(1, lambda x: x[0]<=length/2),\ + (2, lambda x: np.logical_and(x[0]>=length/2, x[1]<=height/2)), + (3, lambda x: np.logical_and(x[0]>=length/2, x[1]>=height/2))] + + cell_tags = make_cell_tags(domain, subdomains) + facet_tags = make_facet_tags(domain, boundaries) + + * Common **material laws**: - .. code-block:: python + .. jupyter-execute:: # V is a dolfinx.fem.function_space # cell_tags is a dolfinx.mesh.MeshTags object @@ -30,7 +61,7 @@ Using the ``elastodynamicsx.pde`` package: tag_mat1 = 1 # suppose this refers to cells associated with material no 1 tags_mat2 = (2, 3) # same for material no 2 mat1 = material((V, cell_tags, tag_mat1), 'isotropic', rho=1, lambda_=2, mu=1) - mat2 = material((V, cell_tags, tag_mat2), 'isotropic', rho=2, lambda_=4, mu=2) + mat2 = material((V, cell_tags, tags_mat2), 'isotropic', rho=2, lambda_=4, mu=2) * linear elasticity: *scalar*, *isotropic*, *cubic*, *hexagonal*, *trigonal*, *tetragonal*, *orthotropic*, *monoclinic*, *triclinic* @@ -43,7 +74,7 @@ Using the ``elastodynamicsx.pde`` package: * Common **boundary conditions** (BCs): - .. code-block:: python + .. jupyter-execute:: # facet_tags is a dolfinx.mesh.MeshTags object @@ -62,21 +93,23 @@ Using the ``elastodynamicsx.pde`` package: * Define **body forces**: - .. code-block:: python + .. jupyter-execute:: from elastodynamicsx.pde import BodyForce amplitude = default_scalar_type([0, 0]) # a dummy load amplitude def shape_x(x): - return np.ones_like(x) # a dummy shape + x1, x2 = 0.2, 0.3 + y1, y2 = 0.4, 0.5 + return (x[0] >= x1) * (x[0] <= x2) * (x[1] >= y1) * (x[1] <= y2) # a dummy shape f_body = fem.Function(V) f_body.interpolate(lambda x: amplitude[:,np.newaxis] * shape_x(x)[np.newaxis,:]) - f1 = BodyForce(V, f_body) # not specifying cell_tags and a specific tag means the entire domain + f1 = BodyForce((V, cell_tags, None), f_body) # None for the entire domain * **Assemble** several materials, BCs and body forces into a *PDE* instance: - .. code-block:: python + .. jupyter-execute:: from elastodynamicsx.pde import PDE @@ -106,7 +139,8 @@ Using the ``elastodynamicsx.solvers`` package: .. tab:: Time domain - .. code-block:: python + .. jupyter-execute:: + :hide-output: # Time integration from elastodynamicsx.solvers import TimeStepper @@ -119,7 +153,7 @@ Using the ``elastodynamicsx.solvers`` package: T_N.value = np.sin(t) * forceVector # Initialize the time stepper: compile forms and assemble the mass matrix - tStepper = TimeStepper.build(V, pde.m, pde.c, pde.k, pde.L, dt, bcs=bcs, scheme='leapfrog') + tStepper = TimeStepper.build(V, pde.m, pde.c, pde.k, pde.L, dt, bcs=pde.bcs, scheme='newmark') # Define the initial values tStepper.set_initial_condition(u0=[0, 0], v0=[0, 0], t0=0) @@ -169,7 +203,8 @@ Using the ``elastodynamicsx.solvers`` package: .. tab:: Eigenmodes - .. code-block:: python + .. jupyter-execute:: + :hide-output: # Normal modes from elastodynamicsx.solvers import EigenmodesSolver @@ -217,7 +252,8 @@ Using the ``elastodynamicsx.solutions`` package: * **Eigenmodes** solutions: -.. code-block:: python +.. jupyter-execute:: + :hide-output: # eps is a elastodynamicsx.solvers.EigenmodesSolver # eps.solve() has already been performed From 4492cfed104f590498887707597c0587a67d622b Mon Sep 17 00:00:00 2001 From: Pierric Mora Date: Fri, 5 Apr 2024 16:14:20 +0200 Subject: [PATCH 04/12] pytest and debug spectral elements --- elastodynamicsx/utils/spectralelements.py | 30 ++++++++--- test/utils/test_spectralelements.py | 66 +++++++++++++++++++++++ 2 files changed, 89 insertions(+), 7 deletions(-) create mode 100644 test/utils/test_spectralelements.py diff --git a/elastodynamicsx/utils/spectralelements.py b/elastodynamicsx/utils/spectralelements.py index 9db3d5a..f32fda3 100644 --- a/elastodynamicsx/utils/spectralelements.py +++ b/elastodynamicsx/utils/spectralelements.py @@ -46,7 +46,11 @@ def GLL_element(cell_type, shape: typing.Optional[typing.Tuple[int, ...]] = None) -> basix.ufl._BasixElement: """Element defined using the Gauss-Lobatto-Legendre points""" cell_type = _suitable_cell_type_format(cell_type) - element = basix.ufl.element(basix.ElementFamily.P, cell_type, degree, basix.LagrangeVariant.gll_warped, shape=shape) + element = basix.ufl.element(basix.ElementFamily.P, + cell_type, + degree, + basix.LagrangeVariant.gll_warped, + shape=shape) return element @@ -55,7 +59,12 @@ def GL_element(cell_type, shape: typing.Optional[typing.Tuple[int, ...]] = None) -> basix.ufl._BasixElement: """(discontinuous) Element defined using the Gauss-Legendre points""" cell_type = _suitable_cell_type_format(cell_type) - element = basix.ufl.element(basix.ElementFamily.P, cell_type, degree, basix.LagrangeVariant.gl_warped, True) + element = basix.ufl.element(basix.ElementFamily.P, + cell_type, + degree, + basix.LagrangeVariant.gl_warped, + discontinuous=True, + shape=shape) return element @@ -64,7 +73,12 @@ def Legendre_element(cell_type, shape: typing.Optional[typing.Tuple[int, ...]] = None) -> basix.ufl._BasixElement: """(discontinuous) Element whose basis functions are the orthonormal Legendre polynomials""" cell_type = _suitable_cell_type_format(cell_type) - element = basix.ufl.element(basix.ElementFamily.P, cell_type, degree, basix.LagrangeVariant.legendre, True) + element = basix.ufl.element(basix.ElementFamily.P, + cell_type, + degree, + basix.LagrangeVariant.legendre, + discontinuous=True, + shape=shape) return element @@ -74,20 +88,22 @@ def Legendre_element(cell_type, def GLL_quadrature(degree: int) -> dict: """Returns the dolfinx quadrature rule for use with GLL elements of the given degree.""" - if degree == 2: - return {"quadrature_rule": "GLL", "quadrature_degree": 3} + if degree <= 2: + return {"quadrature_rule": "GLL", "quadrature_degree": degree} else: return {"quadrature_rule": "GLL", "quadrature_degree": 2 * (degree - 1)} def GL_quadrature(degree: int) -> dict: """Returns the dolfinx quadrature rule for use with GL elements of the given degree.""" - return {"quadrature_rule": "GL", "quadrature_degree": 2 * degree} # 2 * degree + 1 ? + return {"quadrature_rule": "GL", "quadrature_degree": 2 * degree} def Legendre_quadrature(degree: int) -> dict: """Returns the dolfinx quadrature rule for use with Legendre elements of the given degree.""" - return {"quadrature_rule": "GL", "quadrature_degree": 2 * degree} + if degree <= 2: + return {"quadrature_rule": "GL", "quadrature_degree": degree} + return {"quadrature_rule": "GL", "quadrature_degree": 2 * (degree - 1)} # ## ### ### ### ### diff --git a/test/utils/test_spectralelements.py b/test/utils/test_spectralelements.py new file mode 100644 index 0000000..b19eccb --- /dev/null +++ b/test/utils/test_spectralelements.py @@ -0,0 +1,66 @@ +# Copyright (C) 2024 Pierric Mora +# +# This file is part of ElastodynamiCSx +# +# SPDX-License-Identifier: MIT + +import numpy as np + +from mpi4py import MPI +from petsc4py import PETSc +from dolfinx import mesh, fem +import ufl + +from elastodynamicsx.pde import PDE, material, PDECONFIG +from elastodynamicsx.utils import spectral_element, spectral_quadrature + + +def test(): + nbelts = 3 + domain = mesh.create_unit_square(MPI.COMM_WORLD, nbelts, nbelts, + cell_type=mesh.CellType.quadrilateral) + + degrees = range(1, 10) + names = ('GLL', 'GL', 'LEGENDRE') + for degree in degrees: + for name in names: + specFE = spectral_element(name, mesh.CellType.quadrilateral, degree) + specmd = spectral_quadrature(name, degree) + V = fem.FunctionSpace(domain, specFE) + + ####### + # Compile mass matrices using pure fenicsx code + u, v = ufl.TrialFunction(V), ufl.TestFunction(V) + + # here we specify the quadrature metadata -> the matrix becomes diagonal + a1_diag = fem.form(ufl.inner(u, v) * ufl.dx(metadata=specmd)) + M1_diag = fem.petsc.assemble_matrix(a1_diag) + M1_diag.assemble() + + ####### + # Compile mass matrices using elastodynamicsx.pde classes + # elastodynamicsx offers two ways to specify the quadrature metadata + + # 1. by passing the metadata as kwargs to each material / boundary condition / bodyforce + + mat2 = material(V, 'scalar', 1, 1, metadata=specmd) + pde2_diag = PDE(V, [mat2]) + M2_diag = pde2_diag.M() + + # 2. by passing the metadata to the PDE.default_metadata, which is being used by + # all materials / boundary conditions / bodyforces at __init__ call + + PDECONFIG.default_metadata = specmd # default for all forms in the pde package (materials, BCs, ...) + + mat3 = material(V, 'scalar', 1, 1) + pde3_diag = PDE(V, [mat3]) + M3_diag = pde3_diag.M() + + for M_ in (M1_diag, M2_diag, M3_diag): + d_ = M_.getDiagonal() + MforceDiag = PETSc.Mat() + MforceDiag.createAIJ(M_.size) + MforceDiag.assemble() + MforceDiag.setDiagonal(d_) + assert np.isclose((M_ - MforceDiag).norm(), 0.), \ + f"Matrix should be diagonal; Something went wrong for elts:{name}, deg:{degree}" From bd22151c9cbb0705f354a17881e5f449096b098b Mon Sep 17 00:00:00 2001 From: Pierric Mora Date: Mon, 8 Apr 2024 11:05:14 +0200 Subject: [PATCH 05/12] dolfinx v0.7.3 --- Dockerfile.lab | 2 +- Dockerfile.shell | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Dockerfile.lab b/Dockerfile.lab index b1494e4..1b7d29f 100644 --- a/Dockerfile.lab +++ b/Dockerfile.lab @@ -1,4 +1,4 @@ -FROM dolfinx/lab:v0.7.2 +FROM dolfinx/lab:v0.7.3 LABEL maintainer="Pierric Mora " LABEL description=" elastodynamics with FEniCSx/DOLFINx " diff --git a/Dockerfile.shell b/Dockerfile.shell index 7807619..d311116 100644 --- a/Dockerfile.shell +++ b/Dockerfile.shell @@ -1,4 +1,4 @@ -FROM dolfinx/dolfinx:v0.7.2 +FROM dolfinx/dolfinx:v0.7.3 LABEL maintainer="Pierric Mora " LABEL description=" elastodynamics with FEniCSx/DOLFINx " From 7596ace2aff778c1df383fa824e0c2524a3ac94e Mon Sep 17 00:00:00 2001 From: Pierric Mora Date: Fri, 12 Apr 2024 15:11:16 +0200 Subject: [PATCH 06/12] dynamic version info for dolfinx_mpc --- docsrc/gettingstarted/installation.rst | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/docsrc/gettingstarted/installation.rst b/docsrc/gettingstarted/installation.rst index 50268b0..c0bdddd 100644 --- a/docsrc/gettingstarted/installation.rst +++ b/docsrc/gettingstarted/installation.rst @@ -89,7 +89,13 @@ Dependencies import dolfinx print(f"DOLFINx version: {dolfinx.__version__}") -* `DOLFINx-MPC `_ v0.7.0.post1. This dependency is optional (periodic BCs). +* `DOLFINx-MPC `_. This dependency is optional (periodic BCs). + +.. jupyter-execute:: + :hide-code: + + from importlib.metadata import version + print(f"DOLFINx-MPC version: {version('dolfinx_mpc')}") * ``numpy`` From 08f21e1f1959ce7a0af9ce3e25943c8f60fcabb4 Mon Sep 17 00:00:00 2001 From: Pierric Mora Date: Mon, 22 Apr 2024 17:16:13 +0200 Subject: [PATCH 07/12] ignore anim outputs in demo/ --- .gitignore | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index 7d326ef..c351d4d 100644 --- a/.gitignore +++ b/.gitignore @@ -6,7 +6,8 @@ # *.JPG *.npz # *.npy -# *.mp4 +demo/*.gif +demo/*.mp4 # *.html # *.sqlite # *.hdf5 From 34e4a33ae93e5d3ac863646d98a408fbd6c48ab3 Mon Sep 17 00:00:00 2001 From: Pierric Mora Date: Mon, 22 Apr 2024 17:51:49 +0200 Subject: [PATCH 08/12] export anims with live_plotter --- elastodynamicsx/plot/plot.py | 68 +++++++++++++++++-- .../frequencydomain/frequencydomainsolver.py | 9 ++- .../solvers/timedomain/timesteppers.py | 7 +- 3 files changed, 78 insertions(+), 6 deletions(-) diff --git a/elastodynamicsx/plot/plot.py b/elastodynamicsx/plot/plot.py index 70b094a..43e0491 100644 --- a/elastodynamicsx/plot/plot.py +++ b/elastodynamicsx/plot/plot.py @@ -278,6 +278,32 @@ def update_scalars(self, *all_scalars, **kwargs): if kwargs.get('render', True): self.render() + def live_plotter_start(self): + is_recording = hasattr(self, 'mwriter') + if is_recording: + self.notebook = False + self.show(interactive_update=True) + + def live_plotter_stop(self): + is_recording = hasattr(self, 'mwriter') + if is_recording: + self.close() + fname = self.mwriter.request.filename + try: + import IPython.display + if self.mwriter.request.extension.lower() == '.gif': + im = IPython.display.Image(open(fname, 'rb').read()) + IPython.display.display(im) # display .gif in notebook + elif self.mwriter.request.extension.lower() == '.mp4': + from elastodynamicsx import _DOCS_CFG + vi = IPython.display.Video(fname) + if _DOCS_CFG: + vi.embed = True + vi.html_attributes = "controls loop autoplay" + IPython.display.display(vi) # display .mp4 in notebook + except ModuleNotFoundError: + pass + def live_plotter_update_function(self, i: int, vec: PETSc.Vec) -> None: # type: ignore[name-defined] if not isinstance(vec, PETSc.Vec): # type: ignore[attr-defined] if issubclass(type(vec), fem.function.Function): @@ -289,9 +315,13 @@ def live_plotter_update_function(self, i: int, vec: PETSc.Vec) -> None: # type: except: # noqa raise TypeError if (self._refresh_step > 0) and (i % self._refresh_step == 0): + is_recording = hasattr(self, 'mwriter') with vec.localForm() as loc_v: # Necessary for correct handling of ghosts in parallel - self.update_scalars(loc_v.array) - time.sleep(self._tsleep) + self.update_scalars(loc_v.array, render=not is_recording) + if is_recording: + self.write_frame() # This triggers a render + else: + time.sleep(self._tsleep) def add_time_browser(self, update_fields_function: Callable, @@ -436,6 +466,32 @@ def update_vectors(self, *all_vectors, render: bool = True) -> None: if render: self.render() + def live_plotter_start(self): + is_recording = hasattr(self, 'mwriter') + if is_recording: + self.notebook = False + self.show(interactive_update=True) + + def live_plotter_stop(self): + is_recording = hasattr(self, 'mwriter') + if is_recording: + self.close() + fname = self.mwriter.request.filename + try: + import IPython.display + if self.mwriter.request.extension.lower() == '.gif': + im = IPython.display.Image(open(fname, 'rb').read()) + IPython.display.display(im) # display .gif in notebook + elif self.mwriter.request.extension.lower() == '.mp4': + from elastodynamicsx import _DOCS_CFG + vi = IPython.display.Video(fname) + if _DOCS_CFG: + vi.embed = True + vi.html_attributes = "controls loop autoplay" + IPython.display.display(vi) # display .mp4 in notebook + except ModuleNotFoundError: + pass + def live_plotter_update_function(self, i: int, vec: PETSc.Vec) -> None: # type: ignore[name-defined] if not isinstance(vec, PETSc.Vec): # type: ignore[attr-defined] if issubclass(type(vec), fem.function.Function): @@ -446,9 +502,13 @@ def live_plotter_update_function(self, i: int, vec: PETSc.Vec) -> None: # type: except: # noqa raise TypeError if (self._refresh_step > 0) and (i % self._refresh_step == 0): + is_recording = hasattr(self, 'mwriter') with vec.localForm() as loc_v: # Necessary for correct handling of ghosts in parallel - self.update_vectors(loc_v.array) - time.sleep(self._tsleep) + self.update_vectors(loc_v.array, render=not is_recording) + if is_recording: + self.write_frame() # This triggers a render + else: + time.sleep(self._tsleep) def add_time_browser(self, update_fields_function: Callable, diff --git a/elastodynamicsx/solvers/frequencydomain/frequencydomainsolver.py b/elastodynamicsx/solvers/frequencydomain/frequencydomainsolver.py index 0c010bd..ba7f3e3 100644 --- a/elastodynamicsx/solvers/frequencydomain/frequencydomainsolver.py +++ b/elastodynamicsx/solvers/frequencydomain/frequencydomainsolver.py @@ -169,6 +169,8 @@ def _solve_multiple_omegas(self, omegas: np.ndarray, # Loop on values in omegas -> _solve_single_omega + call_on_leave = kwargs.get('call_on_leave', []) + live_plt = kwargs.get('live_plotter', None) if not (live_plt is None): @@ -176,14 +178,19 @@ def _solve_multiple_omegas(self, omegas: np.ndarray, # from elastodynamicsx.plot import live_plotter # live_plt = live_plotter(self.u, live_plt.pop('refresh_step', 1), **live_plt) raise NotImplementedError + callbacks.append(live_plt.live_plotter_update_function) - live_plt.show(interactive_update=True) + call_on_leave.append(live_plt.live_plotter_stop) + live_plt.live_plotter_start() for i in tqdm(range(len(omegas))): self._solve_single_omega(omegas[i], out) for callback in callbacks: callback(i, out) # <- store solution, plot, print, ... + for function in call_on_leave: + function() + return out @property diff --git a/elastodynamicsx/solvers/timedomain/timesteppers.py b/elastodynamicsx/solvers/timedomain/timesteppers.py index dfa19e0..37df6ba 100644 --- a/elastodynamicsx/solvers/timedomain/timesteppers.py +++ b/elastodynamicsx/solvers/timedomain/timesteppers.py @@ -280,6 +280,7 @@ def solve(self, num_steps, **kwargs): callfirsts = kwargs.get('callfirsts', []) callbacks = kwargs.get('callbacks', []) + call_on_leave = kwargs.get('call_on_leave', []) live_plt = kwargs.get('live_plotter', None) if not (live_plt is None): @@ -289,7 +290,8 @@ def solve(self, num_steps, **kwargs): live_plt = live_plotter(u_init, live_plt.pop('refresh_step', 1), **live_plt) callbacks.append(live_plt.live_plotter_update_function) - live_plt.show(interactive_update=True) + call_on_leave.append(live_plt.live_plotter_stop) + live_plt.live_plotter_start() t0 = self.t self._tscheme.initialStep(t0, callfirsts, callbacks, verbose=verbose) @@ -331,3 +333,6 @@ def solve(self, num_steps, **kwargs): for callback in callbacks: callback(i, self._out) # <- store solution, plot, print, ... + + for function in call_on_leave: + function() From 7d9d14912770b682e4ad6baa014ae8ab5cd2fd26 Mon Sep 17 00:00:00 2001 From: Pierric Mora Date: Mon, 22 Apr 2024 17:52:05 +0200 Subject: [PATCH 09/12] mypy --- test/utils/test_spectralelements.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/utils/test_spectralelements.py b/test/utils/test_spectralelements.py index b19eccb..d627926 100644 --- a/test/utils/test_spectralelements.py +++ b/test/utils/test_spectralelements.py @@ -9,7 +9,7 @@ from mpi4py import MPI from petsc4py import PETSc from dolfinx import mesh, fem -import ufl +import ufl # type: ignore from elastodynamicsx.pde import PDE, material, PDECONFIG from elastodynamicsx.utils import spectral_element, spectral_quadrature From f05f3a439811bec74101d1730f243850dd6d9895 Mon Sep 17 00:00:00 2001 From: Pierric Mora Date: Tue, 23 Apr 2024 00:16:10 +0200 Subject: [PATCH 10/12] auto record anim on doc build --- elastodynamicsx/__init__.py | 4 ++-- elastodynamicsx/plot/plot.py | 17 +++++++++++++++-- 2 files changed, 17 insertions(+), 4 deletions(-) diff --git a/elastodynamicsx/__init__.py b/elastodynamicsx/__init__.py index 9a8f91a..d36a714 100644 --- a/elastodynamicsx/__init__.py +++ b/elastodynamicsx/__init__.py @@ -8,8 +8,6 @@ from ._version import __version__ # noqa -from elastodynamicsx import pde, plot, solutions, solvers, utils - import os # A simple flag to switch the jupyter_backend when building the docs _DOCS_CFG = os.environ.get("ELASTODYNAMICSX_DOCS_CFG", "FALSE").upper() == "TRUE" @@ -18,4 +16,6 @@ import pyvista pyvista.set_jupyter_backend('static') +from elastodynamicsx import pde, plot, solutions, solvers, utils # noqa + __all__ = ["pde", "plot", "solutions", "solvers", "utils"] diff --git a/elastodynamicsx/plot/plot.py b/elastodynamicsx/plot/plot.py index 43e0491..03fa3f6 100644 --- a/elastodynamicsx/plot/plot.py +++ b/elastodynamicsx/plot/plot.py @@ -17,6 +17,8 @@ from dolfinx import plot, fem from dolfinx.mesh import Mesh, MeshTags +from elastodynamicsx import _DOCS_CFG + # ## ------------------------------------------------------------------------- ## # # ## --- preliminary: auto-configure pyvista backend for jupyter notebooks --- ## # @@ -295,7 +297,6 @@ def live_plotter_stop(self): im = IPython.display.Image(open(fname, 'rb').read()) IPython.display.display(im) # display .gif in notebook elif self.mwriter.request.extension.lower() == '.mp4': - from elastodynamicsx import _DOCS_CFG vi = IPython.display.Video(fname) if _DOCS_CFG: vi.embed = True @@ -466,7 +467,20 @@ def update_vectors(self, *all_vectors, render: bool = True) -> None: if render: self.render() + def _auto_record(self, fname=None, *args): + is_recording = hasattr(self, 'mwriter') + if is_recording: + return + if fname is None: + fname = 'tmp.mp4' + if fname.endswith('.gif'): + self.open_gif(fname, *args) + elif fname.endswith('.mp4'): + self.open_movie(fname, *args) + def live_plotter_start(self): + if _DOCS_CFG: + self._auto_record() is_recording = hasattr(self, 'mwriter') if is_recording: self.notebook = False @@ -483,7 +497,6 @@ def live_plotter_stop(self): im = IPython.display.Image(open(fname, 'rb').read()) IPython.display.display(im) # display .gif in notebook elif self.mwriter.request.extension.lower() == '.mp4': - from elastodynamicsx import _DOCS_CFG vi = IPython.display.Video(fname) if _DOCS_CFG: vi.embed = True From aaf7c0dc705344231f524fecfba00f3405aca188 Mon Sep 17 00:00:00 2001 From: Pierric Mora Date: Tue, 23 Apr 2024 12:00:29 +0200 Subject: [PATCH 11/12] update demos with anims --- demo/eigs_3D_AluminumCube.py | 18 ++-- demo/eigs_3D_ElasticBeam.py | 30 +++---- demo/freq_2D-PSV_FullSpace.py | 52 +++++++----- demo/freq_2D-SH_FullSpace.py | 51 +++++------ demo/tdsdyn_3D_ElasticBeam.py | 64 +++++++------- demo/weq_2D-PSV_FullSpace.py | 55 +++++++----- ...lfSpace_Lamb_KomatitschVilotte_BSSA1998.py | 49 ++++++----- demo/weq_2D-SH_FullSpace.py | 70 +++++++++------- demo/weqnl_1D-PSV_Murnaghan_Pwave.py | 84 +++++++++++-------- 9 files changed, 263 insertions(+), 210 deletions(-) diff --git a/demo/eigs_3D_AluminumCube.py b/demo/eigs_3D_AluminumCube.py index 048953f..e929c27 100644 --- a/demo/eigs_3D_AluminumCube.py +++ b/demo/eigs_3D_AluminumCube.py @@ -23,7 +23,6 @@ from dolfinx import mesh, fem, default_scalar_type from mpi4py import MPI -from petsc4py import PETSc from elastodynamicsx.pde import material, PDE from elastodynamicsx.solvers import EigenmodesSolver @@ -83,7 +82,7 @@ eigenfreqs = eps.getEigenfrequencies() # eigenmodes = eps.getEigenmodes() -eps.plot(V, slice(6,6+9), wireframe=True, factor=30) # Avoid the first 6 rigid body modes +eps.plot(V, slice(6, 6+9), wireframe=True, factor=30) # Avoid the first 6 rigid body modes # - @@ -102,14 +101,15 @@ 310.109, 316.197, 317.392, 326.462, 329.034, 332.441, 333.364, 336.65, 337.359, 338.276]) -freqs_OgiEtAl_calc = np.array([116.32 , 143.186, 158.44 , 166.113, 169.338, 178.36 , 184.57 , 185.078, \ - 190.206, 197.692, 201.462, 207.096, 211 , 215.613, 223.219, 230.804, \ - 233.329, 234.758, 250.777, 251.038, 252.303, 256.849, 258.064, 258.874, \ - 259.203, 267.746, 276.736, 279.144, 282.773, 293.016, 304.593, 305.316, \ - 309.591, 315.775, 317.931, 326.556, 329.369, 332.732, 332.271, 336.218, \ +freqs_OgiEtAl_calc = np.array([116.32 , 143.186, 158.44 , 166.113, 169.338, 178.36 , 184.57 , 185.078, + 190.206, 197.692, 201.462, 207.096, 211 , 215.613, 223.219, 230.804, + 233.329, 234.758, 250.777, 251.038, 252.303, 256.849, 258.064, 258.874, + 259.203, 267.746, 276.736, 279.144, 282.773, 293.016, 304.593, 305.316, + 309.591, 315.775, 317.931, 326.556, 329.369, 332.732, 332.271, 336.218, 337.511, 337.71]) print('Eigenfrequencies: comparison with litterature values') print(' FE \tOgi et al, calc.\t Ogi et al, exp. \t(kHz)') -for fFE, fOgi_calc, fOgi_exp in zip(eigenfreqs[6:]*1e3, freqs_OgiEtAl_calc, freqs_OgiEtAl_exp): # *1e3 to convert MHz into kHz - print(str(round(fFE, 3)) +"\t "+ str(round(fOgi_calc, 3)) +"\t\t "+ str(round(fOgi_exp, 3))) +for fFE, fOgi_calc, fOgi_exp in zip(eigenfreqs[6:] * 1e3, freqs_OgiEtAl_calc, freqs_OgiEtAl_exp): + # *1e3 to convert MHz into kHz + print(f"{fFE:.3f}\t {fOgi_calc:.3f}\t\t {fOgi_exp:.3f}") diff --git a/demo/eigs_3D_ElasticBeam.py b/demo/eigs_3D_ElasticBeam.py index b997e0f..7874eef 100644 --- a/demo/eigs_3D_ElasticBeam.py +++ b/demo/eigs_3D_ElasticBeam.py @@ -30,15 +30,10 @@ from dolfinx import mesh, fem, default_scalar_type from mpi4py import MPI -from petsc4py import PETSc from elastodynamicsx.pde import material, boundarycondition, PDE from elastodynamicsx.solvers import EigenmodesSolver from elastodynamicsx.utils import make_facet_tags - -#import pyvista -#pyvista.start_xvfb() -#pyvista.set_jupyter_backend("static") # - # ### FE domain @@ -48,8 +43,8 @@ # Nb of elts. Nx = 20 -Ny = int(B_/L_ * Nx) + 1 -Nz = int(H_/L_ * Nx) + 1 +Ny = int(B_ / L_ * Nx) + 1 +Nz = int(H_ / L_ * Nx) + 1 extent = [[0., 0., 0.], [L_, B_, H_]] @@ -61,11 +56,11 @@ # define some tags tag_left, tag_top, tag_right, tag_bottom, tag_back, tag_front = 1, 2, 3, 4, 5, 6 -boundaries = [(tag_left , lambda x: np.isclose(x[0], 0 )),\ - (tag_right , lambda x: np.isclose(x[0], L_)),\ - (tag_bottom, lambda x: np.isclose(x[1], 0 )),\ - (tag_top , lambda x: np.isclose(x[1], B_)),\ - (tag_back , lambda x: np.isclose(x[2], 0 )),\ +boundaries = [(tag_left , lambda x: np.isclose(x[0], 0.)), + (tag_right , lambda x: np.isclose(x[0], L_)), + (tag_bottom, lambda x: np.isclose(x[1], 0.)), + (tag_top , lambda x: np.isclose(x[1], B_)), + (tag_back , lambda x: np.isclose(x[2], 0.)), (tag_front , lambda x: np.isclose(x[2], H_))] facet_tags = make_facet_tags(domain, boundaries) @@ -88,7 +83,7 @@ # Scaling scaleRHO = 1e6 # a scaling factor to avoid blowing the solver -scaleFREQ= np.sqrt(scaleRHO) # the frequencies must be scaled accordingly +scaleFREQ = np.sqrt(scaleRHO) # the frequencies must be scaled accordingly rho *= scaleRHO # Convert Young & Poisson to Lamé's constants @@ -134,12 +129,13 @@ # Exact solution computation from scipy.optimize import root from math import cos, cosh -falpha = lambda x: cos(x)*cosh(x)+1 -alpha = lambda n: root(falpha, (2*n+1)*np.pi/2)['x'][0] +falpha = lambda x: cos(x) * cosh(x) + 1 +alpha = lambda n: root(falpha, (2 * n + 1) * np.pi / 2)['x'][0] nev = eigenfreqs.size -I_bend = H_*B_**3/12*(np.arange(nev)%2==0) + B_*H_**3/12*(np.arange(nev)%2==1) -freq_beam = np.array([alpha(i//2) for i in range(nev)])**2 *np.sqrt(E*I_bend/(rho.value*B_*H_*L_**4))/2/np.pi +I_bend = H_ * B_**3 / 12 * (np.arange(nev) % 2 == 0) + B_ * H_**3 / 12 * (np.arange(nev) % 2 == 1) +freq_beam = np.array([alpha(i // 2) for i in range(nev)])**2 \ + * np.sqrt(E * I_bend / (rho.value * B_ * H_ * L_**4)) / 2 / np.pi print('Eigenfrequencies: comparison with beam theory\n') print('mode || FE (Hz)\t|| Beam theory (Hz)\t|| Difference (%)') diff --git a/demo/freq_2D-PSV_FullSpace.py b/demo/freq_2D-PSV_FullSpace.py index 6ac301c..be6e553 100644 --- a/demo/freq_2D-PSV_FullSpace.py +++ b/demo/freq_2D-PSV_FullSpace.py @@ -26,12 +26,11 @@ from dolfinx import mesh, fem, default_scalar_type import ufl from mpi4py import MPI -from petsc4py import PETSc from elastodynamicsx.pde import material, BodyForce, boundarycondition, PDE from elastodynamicsx.solvers import FrequencyDomainSolver from elastodynamicsx.plot import plotter, live_plotter -from elastodynamicsx.utils import make_facet_tags, make_cell_tags, ParallelEvaluator +from elastodynamicsx.utils import make_facet_tags, ParallelEvaluator from analyticalsolutions import u_2D_PSV_xw, int_Fraunhofer_2D @@ -44,7 +43,7 @@ # + degElement = 1 length, height = 10, 10 -Nx, Ny = 100//degElement, 100//degElement +Nx, Ny = 100 // degElement, 100 // degElement # create the mesh extent = [[0., 0.], [length, height]] @@ -55,9 +54,9 @@ tag_left, tag_top, tag_right, tag_bottom = 1, 2, 3, 4 all_tags = (tag_left, tag_top, tag_right, tag_bottom) -boundaries = [(tag_left , lambda x: np.isclose(x[0], 0 )),\ - (tag_right , lambda x: np.isclose(x[0], length)),\ - (tag_bottom, lambda x: np.isclose(x[1], 0 )),\ +boundaries = [(tag_left , lambda x: np.isclose(x[0], 0 )), + (tag_right , lambda x: np.isclose(x[0], length)), + (tag_bottom, lambda x: np.isclose(x[1], 0 )), (tag_top , lambda x: np.isclose(x[1], height))] # define some tags @@ -84,7 +83,7 @@ # + Z_N, Z_T = mat.Z_N, mat.Z_T # P and S mechanical impedances -bc = boundarycondition((V, facet_tags, all_tags), 'Dashpot', Z_N, Z_T) +bc = boundarycondition((V, facet_tags, all_tags), 'Dashpot', Z_N, Z_T) bcs = [bc] # - @@ -96,10 +95,10 @@ # + F0 = fem.Constant(domain, default_scalar_type([1, 0])) # amplitude R0 = 0.1 # radius -X0 = np.array([length/2, height/2, 0]) # center +X0 = np.array([length / 2, height / 2, 0]) # center -x = ufl.SpatialCoordinate(domain) -gaussianBF = F0 * ufl.exp(-((x[0]-X0[0])**2+(x[1]-X0[1])**2)/2/R0**2) / (2*np.pi*R0**2) +x = ufl.SpatialCoordinate(domain) +gaussianBF = F0 * ufl.exp(-((x[0] - X0[0])**2 + (x[1] - X0[1])**2) / 2 / R0**2) / (2 * np.pi * R0**2) bf = BodyForce(V, gaussianBF) # - @@ -156,9 +155,10 @@ from scipy.spatial.transform import Rotation as R theta = np.radians(35) pts = np.linspace(0, length / 2, endpoint=False)[1:] -points_out = X0[:,np.newaxis] + R.from_rotvec([0, 0, theta]).as_matrix() @ np.array([pts, - np.zeros_like(pts), - np.zeros_like(pts)]) +points_out = X0[:, np.newaxis] + \ + R.from_rotvec([0, 0, theta]).as_matrix() @ np.array([pts, + np.zeros_like(pts), + np.zeros_like(pts)]) # Declare a convenience ParallelEvaluator paraEval = ParallelEvaluator(domain, points_out) @@ -167,19 +167,25 @@ u_at_pts_local = np.zeros((paraEval.nb_points_local, V.num_sub_spaces, omegas.size), - dtype=default_scalar_type) # <- output stored here + dtype=default_scalar_type) # <- output stored here + # Callback function: post process solution def cbck_storeAtPoints(i, out): if paraEval.nb_points_local > 0: - u_at_pts_local[:,:,i] = u_res.eval(paraEval.points_local, paraEval.cells_local) + u_at_pts_local[:, :, i] = u_res.eval(paraEval.points_local, paraEval.cells_local) + # Live plotting -enable_plot = True -if domain.comm.rank == 0 and enable_plot: - p = live_plotter(u_res, clim=0.25 * np.linalg.norm(mu.value * F0.value) * np.array([0, 1])) +if domain.comm.rank == 0: + p = live_plotter(u_res, + show_edges=False, + clim=0.25 * np.linalg.norm(mu.value * F0.value) * np.array([0, 1])) if paraEval.nb_points_local > 0: p.add_points(paraEval.points_local) # add points to live_plotter + if p.off_screen: + p.window_size = [640, 480] + p.open_movie('freq_2D-PSV_FullSpace.mp4', framerate=1) else: p = None @@ -196,25 +202,25 @@ def cbck_storeAtPoints(i, out): u_at_pts = paraEval.gather(u_at_pts_local, root=0) if domain.comm.rank == 0: - ### -> Exact solution, At few points + # -> Exact solution, At few points x = points_out.T # account for the size of the source in the analytical formula fn_kdomain_finite_size = int_Fraunhofer_2D['gaussian'](R0) - u_at_pts_anal = u_2D_PSV_xw(x-X0[np.newaxis,:], omegas, F0.value, rho.value, + u_at_pts_anal = u_2D_PSV_xw(x - X0[np.newaxis, :], omegas, F0.value, rho.value, lambda_.value, mu.value, fn_kdomain_finite_size) # fn = np.real icomp = 0 - fig, ax = plt.subplots(len(omegas),1) + fig, ax = plt.subplots(len(omegas), 1) fig.suptitle(r'u at few points, $\theta$=' + str(int(round(np.degrees(theta), 0))) + r'$^{\circ}$') - r = np.linalg.norm(x - X0[np.newaxis,:], axis=1) + r = np.linalg.norm(x - X0[np.newaxis, :], axis=1) for i in range(len(omegas)): ax[i].text(0.15, 0.95, r'$\omega$=' + str(round(omegas[i], 2)), ha='left', va='top', transform=ax[i].transAxes) - ax[i].plot(r, fn(u_at_pts[:, icomp, i]), ls='-' , label='FEM') + ax[i].plot(r, fn(u_at_pts[:, icomp, i]), ls='-', label='FEM') ax[i].plot(r, fn(u_at_pts_anal[:, icomp, i]), ls='--', label='analytical') ax[0].legend() ax[-1].set_xlabel('Distance to source') diff --git a/demo/freq_2D-SH_FullSpace.py b/demo/freq_2D-SH_FullSpace.py index 24b2936..4affc6b 100644 --- a/demo/freq_2D-SH_FullSpace.py +++ b/demo/freq_2D-SH_FullSpace.py @@ -29,13 +29,11 @@ from dolfinx import mesh, fem, default_scalar_type import ufl from mpi4py import MPI -from petsc4py import PETSc from elastodynamicsx.pde import material, BodyForce, boundarycondition, PDE from elastodynamicsx.solvers import FrequencyDomainSolver from elastodynamicsx.plot import plotter, live_plotter -from elastodynamicsx.utils import make_facet_tags, make_cell_tags, ParallelEvaluator -from analyticalsolutions import green_2D_SH_rw, int_Fraunhofer_2D +from elastodynamicsx.utils import make_facet_tags, ParallelEvaluator assert np.issubdtype(default_scalar_type, np.complexfloating), \ "Demo should only be executed with DOLFINx complex mode" @@ -47,20 +45,20 @@ # + degElement = 1 length, height = 10, 10 -Nx, Ny = 100//degElement, 100//degElement +Nx, Ny = 100 // degElement, 100 // degElement # create the mesh extent = [[0., 0.], [length, height]] domain = mesh.create_rectangle(MPI.COMM_WORLD, extent, [Nx, Ny], mesh.CellType.triangle) # create the function space -V = fem.FunctionSpace(domain, ("Lagrange", degElement)) +V = fem.FunctionSpace(domain, ("Lagrange", degElement)) tag_left, tag_top, tag_right, tag_bottom = 1, 2, 3, 4 all_tags = (tag_left, tag_top, tag_right, tag_bottom) -boundaries = [(tag_left , lambda x: np.isclose(x[0], 0 )),\ - (tag_right , lambda x: np.isclose(x[0], length)),\ - (tag_bottom, lambda x: np.isclose(x[1], 0 )),\ +boundaries = [(tag_left , lambda x: np.isclose(x[0], 0 )), + (tag_right , lambda x: np.isclose(x[0], length)), + (tag_bottom, lambda x: np.isclose(x[1], 0 )), (tag_top , lambda x: np.isclose(x[1], height))] # define some tags @@ -162,30 +160,35 @@ from scipy.spatial.transform import Rotation as R theta = np.radians(35) pts = np.linspace(0, length / 2, endpoint=False)[1:] -points_out = X0[:,np.newaxis] + R.from_rotvec([0, 0, theta]).as_matrix() @ np.array([pts, - np.zeros_like(pts), - np.zeros_like(pts)]) +points_out = X0[:, np.newaxis] + \ + R.from_rotvec([0, 0, theta]).as_matrix() @ np.array([pts, + np.zeros_like(pts), + np.zeros_like(pts)]) # Declare a convenience ParallelEvaluator paraEval = ParallelEvaluator(domain, points_out) # Declare data (local) -u_at_pts_local = np.zeros((paraEval.nb_points_local, - 1, - omegas.size), - dtype=default_scalar_type) # <- output stored here +u_at_pts_local = np.zeros((paraEval.nb_points_local, 1, omegas.size), + dtype=default_scalar_type) # <- output stored here + # Callback function: post process solution def cbck_storeAtPoints(i, out): if paraEval.nb_points_local > 0: - u_at_pts_local[:,:,i] = u_res.eval(paraEval.points_local, paraEval.cells_local) + u_at_pts_local[:, :, i] = u_res.eval(paraEval.points_local, paraEval.cells_local) + # Live plotting -enable_plot = True -if domain.comm.rank == 0 and enable_plot: - p = live_plotter(u_res, clim=0.25 * np.linalg.norm(mu.value * F0.value) * np.array([-1, 1])) +if domain.comm.rank == 0: + p = live_plotter(u_res, + show_edges=False, + clim=0.25 * np.linalg.norm(mu.value * F0.value) * np.array([-1, 1])) if paraEval.nb_points_local > 0: p.add_points(paraEval.points_local) # add points to live_plotter + if p.off_screen: + p.window_size = [640, 480] + p.open_movie('freq_2D-SH_FullSpace.mp4', framerate=1) else: p = None @@ -202,27 +205,27 @@ def cbck_storeAtPoints(i, out): u_at_pts = paraEval.gather(u_at_pts_local, root=0) if domain.comm.rank == 0: - ### -> Exact solution, At few points + # -> Exact solution, At few points x = points_out.T - r = np.linalg.norm(x - X0[np.newaxis,:], axis=1) + r = np.linalg.norm(x - X0[np.newaxis, :], axis=1) # account for the size of the source in the analytical formula + from analyticalsolutions import green_2D_SH_rw, int_Fraunhofer_2D fn_kdomain_finite_size = int_Fraunhofer_2D['gaussian'](R0) u_at_pts_anal = green_2D_SH_rw(r, omegas, rho.value, mu.value, fn_kdomain_finite_size) # fn = np.real - fig, ax = plt.subplots(len(omegas),1) + fig, ax = plt.subplots(len(omegas), 1) fig.suptitle(r'u at few points, $\theta$=' + str(int(round(np.degrees(theta), 0))) + r'$^{\circ}$') for i in range(len(omegas)): ax[i].text(0.15, 0.95, r'$\omega$=' + str(round(omegas[i], 2)), ha='left', va='top', transform=ax[i].transAxes) - ax[i].plot(r, fn(u_at_pts[:, 0, i]), ls='-' , label='FEM') + ax[i].plot(r, fn(u_at_pts[:, 0, i]), ls='-', label='FEM') ax[i].plot(r, fn(u_at_pts_anal[:, i]), ls='--', label='analytical') ax[0].legend() ax[-1].set_xlabel('Distance to source') plt.show() # # ------------------ end of Ex 2 ---------------------- - diff --git a/demo/tdsdyn_3D_ElasticBeam.py b/demo/tdsdyn_3D_ElasticBeam.py index fb7d22a..d3893a2 100644 --- a/demo/tdsdyn_3D_ElasticBeam.py +++ b/demo/tdsdyn_3D_ElasticBeam.py @@ -27,11 +27,9 @@ import numpy as np from dolfinx import mesh, fem, default_scalar_type -import ufl from mpi4py import MPI -from petsc4py import PETSc -from elastodynamicsx.pde import material, BodyForce, boundarycondition, PDE, damping +from elastodynamicsx.pde import material, boundarycondition, PDE, damping from elastodynamicsx.solvers import TimeStepper from elastodynamicsx.utils import make_facet_tags, ParallelEvaluator # - @@ -52,11 +50,11 @@ # define some tags tag_left, tag_top, tag_right, tag_bottom, tag_back, tag_front = 1, 2, 3, 4, 5, 6 -boundaries = [(tag_left , lambda x: np.isclose(x[0], 0 )),\ - (tag_right , lambda x: np.isclose(x[0], L_)),\ - (tag_bottom, lambda x: np.isclose(x[1], 0 )),\ - (tag_top , lambda x: np.isclose(x[1], B_)),\ - (tag_back , lambda x: np.isclose(x[2], 0 )),\ +boundaries = [(tag_left , lambda x: np.isclose(x[0], 0.)), + (tag_right , lambda x: np.isclose(x[0], L_)), + (tag_bottom, lambda x: np.isclose(x[1], 0.)), + (tag_top , lambda x: np.isclose(x[1], B_)), + (tag_back , lambda x: np.isclose(x[2], 0.)), (tag_front , lambda x: np.isclose(x[2], H_))] facet_tags = make_facet_tags(domain, boundaries) @@ -68,10 +66,10 @@ # + # Left: clamp -bc_l = boundarycondition((V, facet_tags, tag_left ), 'Clamp') +bc_l = boundarycondition((V, facet_tags, tag_left), 'Clamp') # Right: normal traction (Neumann boundary condition) -T_N = fem.Constant(domain, np.array([0]*3, dtype=default_scalar_type)) +T_N = fem.Constant(domain, np.array([0] * 3, dtype=default_scalar_type)) bc_r = boundarycondition((V, facet_tags, tag_right), 'Neumann', T_N) bcs = [bc_l, bc_r] # list all BCs @@ -82,9 +80,9 @@ # + # -> Time function -p0 = 1. # max amplitude -F_0 = p0 * np.array([0, 0, 1], dtype=default_scalar_type) # source orientation -cutoff_Tc = 4/5 # release time +p_0 = 1. # max amplitude +F_0 = p_0 * np.array([0, 0, 1], dtype=default_scalar_type) # source orientation +cutoff_Tc = 4 / 5 # release time src_t = lambda t: t / cutoff_Tc * (t > 0) * (t <= cutoff_Tc) T_N_function = lambda t: src_t(t) * F_0 @@ -122,7 +120,10 @@ eta_m = fem.Constant(domain, default_scalar_type(eta_m)) eta_k = fem.Constant(domain, default_scalar_type(eta_k)) -material = material(V, 'isotropic', rho, lambda_, mu, damping=damping('Rayleigh', eta_m, eta_k)) +material = material(V, + 'isotropic', + rho, lambda_, mu, + damping=damping('Rayleigh', eta_m, eta_k)) # - # ### Assemble the PDE @@ -141,9 +142,9 @@ # + # Temporal parameters -T = 4 # duration; difference with original example: here t=[0,T-dt] -Nsteps = 50 # number of time steps -dt = T/Nsteps # time increment +T = 4 # duration; difference with original example: here t=[0,T-dt] +Nsteps = 50 # number of time steps +dt = T / Nsteps # time increment # Generalized-alpha method parameters alpha_m = 0.2 @@ -154,7 +155,7 @@ tStepper = TimeStepper.build(V, pde.m, pde.c, pde.k, pde.L, dt, bcs=bcs, **kwargsTScheme) # Set the initial values -tStepper.set_initial_condition(u0=[0,0,0], v0=[0,0,0], t0=0) +tStepper.set_initial_condition(u0=[0, 0, 0], v0=[0, 0, 0], t0=0) # - # ### Define outputs @@ -166,7 +167,7 @@ # + # -> Point evaluation # Define points -points_out = np.array([[L_, B_/2, 0]]).T +points_out = np.array([[L_, B_ / 2, 0]]).T # Declare a convenience ParallelEvaluator paraEval = ParallelEvaluator(domain, points_out) @@ -178,23 +179,25 @@ # -> Energies energies = np.zeros((Nsteps, 4)) -E_damp = 0 # declare; init value +E_damp = 0 # declare; init value u_n = tStepper.timescheme.u # The displacement at time t_n v_n = tStepper.timescheme.v # The velocity at time t_n comm = V.mesh.comm -Energy_elastic = lambda *a: comm.allreduce( fem.assemble_scalar(fem.form( 1/2* pde.k(u_n, u_n) )) , op=MPI.SUM) -Energy_kinetic = lambda *a: comm.allreduce( fem.assemble_scalar(fem.form( 1/2* pde.m(v_n, v_n) )) , op=MPI.SUM) -Energy_damping = lambda *a: dt*comm.allreduce( fem.assemble_scalar(fem.form( pde.c(v_n, v_n) )) , op=MPI.SUM) +Energy_elastic = lambda *a: comm.allreduce(fem.assemble_scalar(fem.form(1/2 * pde.k(u_n, u_n))), op=MPI.SUM) +Energy_kinetic = lambda *a: comm.allreduce(fem.assemble_scalar(fem.form(1/2 * pde.m(v_n, v_n))), op=MPI.SUM) +Energy_damping = lambda *a: dt * comm.allreduce(fem.assemble_scalar(fem.form(pde.c(v_n, v_n))), op=MPI.SUM) # -> live plotting parameters -clim = 0.4 * L_*B_*H_/(E*B_*H_**3/12) * np.amax(F_0)*np.array([0, 1]) -live_plotter = {'refresh_step':1, 'clim':clim} if domain.comm.rank == 0 else None +clim = 0.4 * L_ * B_ * H_ / (E * B_ * H_**3 / 12) * np.amax(F_0) * np.array([0, 1]) +live_plotter = {'refresh_step': 1, 'clim': clim, 'window_size': [640, 480]} if domain.comm.rank == 0 else None + # -> Define callbacks: will be called at the end of each iteration def cbck_storeAtPoints(i, out): if paraEval.nb_points_local > 0: - signals_local[:,:,i+1] = u_n.eval(paraEval.points_local, paraEval.cells_local) + signals_local[:, :, i+1] = u_n.eval(paraEval.points_local, paraEval.cells_local) + def cbck_energies(i, out): global E_damp @@ -202,7 +205,7 @@ def cbck_energies(i, out): E_kin = Energy_kinetic() E_damp+= Energy_damping() E_tot = E_elas + E_kin + E_damp - energies[i+1,:] = np.array([E_elas, E_kin, E_damp, E_tot]) + energies[i+1, :] = np.array([E_elas, E_kin, E_damp, E_tot]) # - @@ -216,8 +219,9 @@ def cbck_energies(i, out): def cfst_updateSources(t): T_N.value = T_N_function(t) + # Run the big time loop! -tStepper.solve(Nsteps-1, +tStepper.solve(Nsteps - 1, callfirsts=[cfst_updateSources], callbacks=[cbck_storeAtPoints, cbck_energies], live_plotter=live_plotter) @@ -233,12 +237,12 @@ def cfst_updateSources(t): # Plot only on rank == 0 if domain.comm.rank == 0: - t = dt*np.arange(energies.shape[0]) + t = dt * np.arange(energies.shape[0]) # Tip displacement u_tip = all_signals[0] plt.figure() - plt.plot(t, u_tip[2,:]) + plt.plot(t, u_tip[2, :]) plt.xlabel('Time') plt.ylabel('Tip displacement') plt.ylim(-0.5, 0.5) diff --git a/demo/weq_2D-PSV_FullSpace.py b/demo/weq_2D-PSV_FullSpace.py index 41517d6..b4c07d6 100644 --- a/demo/weq_2D-PSV_FullSpace.py +++ b/demo/weq_2D-PSV_FullSpace.py @@ -32,7 +32,7 @@ from elastodynamicsx.pde import material, BodyForce, boundarycondition, PDE, PDECONFIG from elastodynamicsx.solvers import TimeStepper from elastodynamicsx.plot import plotter -from elastodynamicsx.utils import spectral_element, spectral_quadrature, make_facet_tags, make_cell_tags, ParallelEvaluator +from elastodynamicsx.utils import spectral_element, spectral_quadrature, make_facet_tags, ParallelEvaluator from analyticalsolutions import u_2D_PSV_xt, int_Fraunhofer_2D # - @@ -52,7 +52,7 @@ # + length, height = 10, 10 -Nx, Ny = 100//degElement, 100//degElement # Nb of elts. +Nx, Ny = 100 // degElement, 100 // degElement # Nb of elts. # create the mesh extent = [[0., 0.], [length, height]] @@ -64,9 +64,9 @@ # define some tags tag_left, tag_top, tag_right, tag_bottom = 1, 2, 3, 4 all_tags = (tag_left, tag_top, tag_right, tag_bottom) -boundaries = [(tag_left , lambda x: np.isclose(x[0], 0 )),\ - (tag_right , lambda x: np.isclose(x[0], length)),\ - (tag_bottom, lambda x: np.isclose(x[1], 0 )),\ +boundaries = [(tag_left , lambda x: np.isclose(x[0], 0 )), + (tag_right , lambda x: np.isclose(x[0], length)), + (tag_bottom, lambda x: np.isclose(x[1], 0 )), (tag_top , lambda x: np.isclose(x[1], height))] facet_tags = make_facet_tags(domain, boundaries) @@ -110,24 +110,30 @@ # Gaussian function nrm = 1 / (2 * np.pi * R0_src**2) # normalize to int[src_x(x) dx]=1 + def src_x(x): # source(x): Gaussian - r = np.linalg.norm(x-X0_src[:,np.newaxis], axis=0) - return nrm * np.exp(-1/2*(r/R0_src)**2, dtype=default_scalar_type) + r = np.linalg.norm(x - X0_src[:, np.newaxis], axis=0) + return nrm * np.exp(-1/2 * (r / R0_src)**2, dtype=default_scalar_type) + # ## -> Time function f0 = 1 # central frequency of the source T0 = 1 / f0 # period d0 = 2 * T0 # duration of source + def src_t(t): # source(t): Sine x Hann window window = np.sin(np.pi * t / d0)**2 * (t < d0) * (t > 0) # Hann window return np.sin(2 * np.pi * f0 * t) * window + # ## -> Space-Time function F_0 = default_scalar_type([1, 0]) # amplitude of the source + def F_body_function(t): # source(x) at a given time - return lambda x: F_0[:,np.newaxis] * src_t(t) * src_x(x)[np.newaxis,:] + return lambda x: F_0[:, np.newaxis] * src_t(t) * src_x(x)[np.newaxis, :] + # Body force 'F_body' F_body = fem.Function(V) # body force @@ -158,10 +164,12 @@ def F_body_function(t): # source(x) at a given time # Time integration # diagonal=True assumes the left hand side operator is indeed diagonal -tStepper = TimeStepper.build(V, pde.m, pde.c, pde.k, pde.L, dt, bcs=bcs, scheme='leapfrog', diagonal=True) +tStepper = TimeStepper.build(V, + pde.m, pde.c, pde.k, pde.L, dt, bcs=bcs, + scheme='leapfrog', diagonal=True) # Set the initial values -tStepper.set_initial_condition(u0=[0,0], v0=[0,0], t0=tstart) +tStepper.set_initial_condition(u0=[0, 0], v0=[0, 0], t0=tstart) # - @@ -170,11 +178,11 @@ def F_body_function(t): # source(x) at a given time # - Live-plot results (only in a terminal; not in a notebook) # + -u_res = tStepper.timescheme.u # The solution +u_res = tStepper.timescheme.u # The solution # -> Extract signals at few points # Define points -points_out = X0_src[:,np.newaxis] + np.array([[1, 1, 0], [2, 2, 0], [3, 3, 0]]).T # shape = (3, nbpts) +points_out = X0_src[:, np.newaxis] + np.array([[1, 1, 0], [2, 2, 0], [3, 3, 0]]).T # shape = (3, nbpts) # Declare a convenience ParallelEvaluator paraEval = ParallelEvaluator(domain, points_out) @@ -184,20 +192,25 @@ def F_body_function(t): # source(x) at a given time V.num_sub_spaces, num_steps)) # <- output stored here + # -> Define callbacks: will be called at the end of each iteration def cbck_storeAtPoints(i, out): if paraEval.nb_points_local > 0: - signals_local[:,:,i+1] = u_res.eval(paraEval.points_local, paraEval.cells_local) + signals_local[:, :, i+1] = u_res.eval(paraEval.points_local, paraEval.cells_local) + # -> enable live plotting enable_plot = True if domain.comm.rank == 0 and enable_plot: - clim = 0.1*np.amax(F_0)*np.array([0, 1]) - kwplot = { 'clim':clim, 'warp_factor':0.5/np.amax(clim) } + clim = 0.1 * np.amax(F_0) * np.array([0, 1]) + kwplot = {'clim': clim, 'warp_factor': 0.5 / np.amax(clim)} p = plotter(u_res, refresh_step=10, **kwplot) if paraEval.nb_points_local > 0: # Add points to live_plotter p.add_points(paraEval.points_local, render_points_as_spheres=True, point_size=12) + if p.off_screen: + p.window_size = [640, 480] + p.open_movie('weq_2D-SH_FullSpace.mp4') else: p = None @@ -213,6 +226,7 @@ def cbck_storeAtPoints(i, out): def cfst_updateSources(t): F_body.interpolate(F_body_function(t)) + # Run the big time loop! tStepper.solve(num_steps - 1, callfirsts=[cfst_updateSources], @@ -235,8 +249,8 @@ def cfst_updateSources(t): # -> Exact solution, At few points x = points_out.T - t = dt*np.arange(num_steps) - signals_exact = u_2D_PSV_xt(x - X0_src[np.newaxis,:], + t = dt * np.arange(num_steps) + signals_exact = u_2D_PSV_xt(x - X0_src[np.newaxis, :], src_t(t), F_0, rho.value, lambda_.value, @@ -244,12 +258,11 @@ def cfst_updateSources(t): dt, fn_kdomain_finite_size) - fig, ax = plt.subplots(1,1) + fig, ax = plt.subplots(1, 1) ax.set_title('Signals at few points') for i in range(len(all_signals)): - ax.plot(t, all_signals[i,0,:], c=f'C{i}', ls='-') # FEM - ax.plot(t, signals_exact[i,0,:], c=f'C{i}', ls='--') # exact + ax.plot(t, all_signals[i, 0, :], c=f'C{i}', ls='-') # FEM + ax.plot(t, signals_exact[i, 0, :], c=f'C{i}', ls='--') # exact ax.set_xlabel('Time') ax.legend(['FEM', 'analytical']) plt.show() - diff --git a/demo/weq_2D-PSV_HalfSpace_Lamb_KomatitschVilotte_BSSA1998.py b/demo/weq_2D-PSV_HalfSpace_Lamb_KomatitschVilotte_BSSA1998.py index 943cff0..2bd0e25 100644 --- a/demo/weq_2D-PSV_HalfSpace_Lamb_KomatitschVilotte_BSSA1998.py +++ b/demo/weq_2D-PSV_HalfSpace_Lamb_KomatitschVilotte_BSSA1998.py @@ -68,6 +68,7 @@ # Create the function space V = fem.FunctionSpace(domain, specFE) + def y_surf(x): """ A convenience function to obtain the 'y' coordinate of a point @@ -88,8 +89,8 @@ def y_surf(x): # + # parameters here... -rho = 2.2 # density -cP, cS = 3.2, 1.8475 # P- and S-wave velocities +rho = 2.2 # density +cP, cS = 3.2, 1.8475 # P- and S-wave velocities # ... end c11, c44 = rho * cP**2, rho * cS**2 @@ -97,7 +98,7 @@ def y_surf(x): mu = fem.Constant(domain, default_scalar_type(c44)) lambda_ = fem.Constant(domain, default_scalar_type(c11 - 2 * c44)) -mat = material(V, 'isotropic', rho, lambda_, mu) +mat = material(V, 'isotropic', rho, lambda_, mu) # - # ### Boundary conditions @@ -126,25 +127,30 @@ def y_surf(x): # Gaussian function nrm = 1 / np.sqrt(2 * np.pi * W0_src**2) # normalize to int[src_x(x) dx]=1 + def src_x(x): # Source(x): Gaussian r = (x[0] - X0_src[0]) / np.cos(np.radians(tilt)) return nrm * np.exp(-1/2 * (r / W0_src)**2, dtype=default_scalar_type) + # ## -> Time function fc = 14.5 # Central frequency sig = np.sqrt(2) / (2 * np.pi * fc) # Gaussian standard deviation -t0 = 4*sig +t0 = 4 * sig + def src_t(t): # Source(t): Ricker - return (1 - ((t-t0)/sig)**2) * np.exp(-0.5*((t-t0)/sig)**2) + return (1 - ((t - t0) / sig)**2) * np.exp(-0.5 * ((t - t0) / sig)**2) + # ## -> Space-Time function p0 = 1. # Max amplitude F_0 = p0 * default_scalar_type([np.sin(np.radians(tilt)), -np.cos(np.radians(tilt))]) # Amplitude of the source + def T_N_function(t): - return lambda x: F_0[:,np.newaxis] * src_t(t) * src_x(x)[np.newaxis,:] # source(x) at a given time + return lambda x: F_0[:, np.newaxis] * src_t(t) * src_x(x)[np.newaxis, :] # source(x) at a given time # - @@ -161,16 +167,18 @@ def T_N_function(t): dt = 0.25e-3 # Time step num_steps = int(6000 * sizefactor) -cmax = ufl.sqrt((lambda_+2*mu)/rho) # max velocity +cmax = ufl.sqrt((lambda_ + 2 * mu) / rho) # max velocity courant_number = TimeStepper.Courant_number(V.mesh, cmax, dt) # Courant number PETSc.Sys.Print(f'CFL condition: Courant number = {courant_number:.2f}') # Time integration # diagonal=True assumes the left hand side operator is indeed diagonal -tStepper = TimeStepper.build(V, pde.m, pde.c, pde.k, pde.L, dt, bcs=bcs, scheme='leapfrog', diagonal=True) +tStepper = TimeStepper.build(V, + pde.m, pde.c, pde.k, pde.L, dt, bcs=bcs, + scheme='leapfrog', diagonal=True) # Set the initial values -tStepper.set_initial_condition(u0=[0,0], v0=[0,0], t0=tstart) +tStepper.set_initial_condition(u0=[0, 0], v0=[0, 0], t0=tstart) # - @@ -197,17 +205,19 @@ def T_N_function(t): V.num_sub_spaces, num_steps)) # <- output stored here + # -> Define callbacks: will be called at the end of each iteration def cbck_storeAtPoints(i, out): if paraEval.nb_points_local > 0: - signals_local[:,:,i+1] = u_res.eval(paraEval.points_local, paraEval.cells_local) + signals_local[:, :, i+1] = u_res.eval(paraEval.points_local, paraEval.cells_local) + # -> enable live plotting enable_plot = True clim = 0.015 * np.linalg.norm(F_0) * np.array([0, 1]) if domain.comm.rank == 0 and enable_plot: - kwplot = {'clim':clim, 'show_edges':False, 'warp_factor':0.05/np.amax(clim) } - p = plotter(u_res, refresh_step=30, **kwplot) + kwplot = {'clim': clim, 'show_edges': False, 'warp_factor': 0.05 / np.amax(clim)} + p = plotter(u_res, refresh_step=30, window_size=[640, 480], **kwplot) if paraEval.nb_points_local > 0: # add points to live_plotter p.add_points(paraEval.points_local, render_points_as_spheres=True, opacity=0.75) @@ -225,6 +235,7 @@ def cbck_storeAtPoints(i, out): def cfst_updateSources(t): T_N.interpolate(T_N_function(t)) + # Run the big time loop! tStepper.solve(num_steps - 1, callfirsts=[cfst_updateSources], @@ -242,8 +253,7 @@ def cfst_updateSources(t): all_signals = paraEval.gather(signals_local, root=0) if domain.comm.rank == 0: - - t = dt*np.arange(num_steps) + t = dt * np.arange(num_steps) # Export as .npz file np.savez('seismogram_weq_2D-PSV_HalfSpace_Lamb_KomatitschVilotte_BSSA1998.npz', @@ -253,14 +263,14 @@ def cfst_updateSources(t): x0 = np.linalg.norm(points_out.T[0] - X0_src) ampl = 4 * dx / np.amax(np.abs(all_signals)) r11, r12 = np.cos(np.radians(tilt)), np.sin(np.radians(tilt)) - # - fig, ax = plt.subplots(1,2) + + fig, ax = plt.subplots(1, 2) ax[0].set_title(r'$u_{tangential}$') ax[1].set_title(r'$u_{normal}$') for i in range(len(all_signals)): - offset = i*dx - x0 - u2plt_t= offset + ampl * ( r11 * all_signals[i,0,:] + r12 * all_signals[i,1,:]) # tangential - u2plt_n= offset + ampl * (-r12 * all_signals[i,0,:] + r11 * all_signals[i,1,:]) # normal + offset = i * dx - x0 + u2plt_t = offset + ampl * ( r11 * all_signals[i, 0, :] + r12 * all_signals[i, 1, :]) # tangential + u2plt_n = offset + ampl * (-r12 * all_signals[i, 0, :] + r11 * all_signals[i, 1, :]) # normal ax[0].plot(u2plt_t, t, c='k') ax[1].plot(u2plt_n, t, c='k') ax[0].fill_betweenx(t, offset, u2plt_t, where=(u2plt_t > offset), color='k') @@ -273,4 +283,3 @@ def cfst_updateSources(t): # TODO: analytical formula - diff --git a/demo/weq_2D-SH_FullSpace.py b/demo/weq_2D-SH_FullSpace.py index 00bc6fd..6741293 100644 --- a/demo/weq_2D-SH_FullSpace.py +++ b/demo/weq_2D-SH_FullSpace.py @@ -33,7 +33,7 @@ from elastodynamicsx.pde import material, BodyForce, boundarycondition, PDE, PDECONFIG from elastodynamicsx.solvers import TimeStepper from elastodynamicsx.plot import plotter -from elastodynamicsx.utils import spectral_element, spectral_quadrature, make_facet_tags, make_cell_tags, ParallelEvaluator +from elastodynamicsx.utils import spectral_element, spectral_quadrature, make_facet_tags, ParallelEvaluator from analyticalsolutions import u_2D_SH_rt, int_Fraunhofer_2D # - @@ -53,21 +53,21 @@ # + length, height = 10, 10 -Nx, Ny = 100//degElement, 100//degElement # Nb of elts. +Nx, Ny = 100 // degElement, 100 // degElement # Nb of elts. # create the mesh extent = [[0., 0.], [length, height]] domain = mesh.create_rectangle(MPI.COMM_WORLD, extent, [Nx, Ny], cell_type) # create the function space -V = fem.FunctionSpace(domain, specFE) +V = fem.FunctionSpace(domain, specFE) # define some tags tag_left, tag_top, tag_right, tag_bottom = 1, 2, 3, 4 all_tags = (tag_left, tag_top, tag_right, tag_bottom) -boundaries = [(tag_left , lambda x: np.isclose(x[0], 0 )),\ - (tag_right , lambda x: np.isclose(x[0], length)),\ - (tag_bottom, lambda x: np.isclose(x[1], 0 )),\ +boundaries = [(tag_left , lambda x: np.isclose(x[0], 0 )), + (tag_right , lambda x: np.isclose(x[0], length)), + (tag_bottom, lambda x: np.isclose(x[1], 0 )), (tag_top , lambda x: np.isclose(x[1], height))] facet_tags = make_facet_tags(domain, boundaries) @@ -106,31 +106,37 @@ # + # ## -> Space function -X0_src = np.array([length/2,height/2,0]) # Center +X0_src = np.array([length / 2, height / 2, 0]) # Center R0_src = 0.1 # Radius # Gaussian function -nrm = 1/(2*np.pi*R0_src**2) # normalize to int[src_x(x) dx]=1 +nrm = 1 / (2 * np.pi * R0_src**2) # normalize to int[src_x(x) dx]=1 + def src_x(x): # source(x): Gaussian - r = np.linalg.norm(x-X0_src[:,np.newaxis], axis=0) - return nrm * np.exp(-1/2*(r/R0_src)**2, dtype=default_scalar_type) + r = np.linalg.norm(x - X0_src[:, np.newaxis], axis=0) + return nrm * np.exp(-1/2 * (r / R0_src)**2, dtype=default_scalar_type) + # ## -> Time function f0 = 1 # central frequency of the source T0 = 1 / f0 # period d0 = 2 * T0 # duration of source + def src_t(t): # source(t): Sine x Hann window - window = np.sin(np.pi*t/d0)**2 * (t0) # Hann window - return np.sin(2*np.pi*f0 * t) * window + window = np.sin(np.pi * t / d0)**2 * (t < d0) * (t > 0) # Hann window + return np.sin(2 * np.pi * f0 * t) * window + # ## -> Space-Time function F_0 = 1 # Amplitude of the source + def F_body_function(t): # source(x) at a given time return lambda x: F_0 * src_t(t) * src_x(x) + # ## Body force 'F_body' F_body = fem.Function(V) # body force gaussianBF = BodyForce(V, F_body) @@ -181,7 +187,7 @@ def F_body_function(t): # source(x) at a given time # -> Extract signals at few points # Define points -points_out = X0_src[:,np.newaxis] + np.array([[1, 0, 0], [2, 0, 0], [3, 0, 0]]).T +points_out = X0_src[:, np.newaxis] + np.array([[1, 0, 0], [2, 0, 0], [3, 0, 0]]).T # Declare a convenience ParallelEvaluator paraEval = ParallelEvaluator(domain, points_out) @@ -196,18 +202,22 @@ def cbck_storeFullField(i, out): if storeAllSteps: all_u[i+1].vector.setArray(out) + def cbck_storeAtPoints(i, out): if paraEval.nb_points_local > 0: - signals_local[:,:,i+1] = u_res.eval(paraEval.points_local, paraEval.cells_local) + signals_local[:, :, i+1] = u_res.eval(paraEval.points_local, paraEval.cells_local) + # -> enable live plotting -enable_plot = True -clim = 0.1*F_0*np.array([-1, 1]) -if domain.comm.rank == 0 and enable_plot: - p = plotter(u_res, refresh_step=10, **{'clim':clim}) +clim = 0.1 * F_0 * np.array([-1, 1]) +if domain.comm.rank == 0: + p = plotter(u_res, refresh_step=10, **{'clim': clim}) if paraEval.nb_points_local > 0: # add points to live_plotter p.add_points(paraEval.points_local, render_points_as_spheres=True, point_size=12) + if p.off_screen: + p.window_size = [640, 480] + p.open_movie('weq_2D-SH_FullSpace.mp4') else: p = None @@ -223,6 +233,7 @@ def cbck_storeAtPoints(i, out): def cfst_updateSources(t): F_body.interpolate(F_body_function(t)) + # Run the big time loop! tStepper.solve(num_steps - 1, callfirsts=[cfst_updateSources], @@ -241,17 +252,17 @@ def cfst_updateSources(t): # Account for the size of the source in the analytical formula fn_kdomain_finite_size = int_Fraunhofer_2D['gaussian'](R0_src) - ### -> Exact solution, Full field + # -> Exact solution, Full field x = u_res.function_space.tabulate_dof_coordinates() - r = np.linalg.norm(x - X0_src[np.newaxis,:], axis=1) - t = dt*np.arange(num_steps) + r = np.linalg.norm(x - X0_src[np.newaxis, :], axis=1) + t = dt * np.arange(num_steps) all_u_n_exact = u_2D_SH_rt(r, src_t(t), rho.value, mu.value, dt, fn_kdomain_finite_size) def update_fields_function(i): - return (all_u[i].x.array, all_u_n_exact[:,i], all_u[i].x.array-all_u_n_exact[:,i]) + return (all_u[i].x.array, all_u_n_exact[:, i], all_u[i].x.array - all_u_n_exact[:, i]) # Initializes with empty fem.Function(V) to have different valid pointers - p = plotter(fem.Function(V), fem.Function(V), fem.Function(V), labels=('FE', 'Exact', 'Diff.'), clim=clim) + p = plotter(fem.Function(V), fem.Function(V), fem.Function(V), labels=('FE', 'Exact', 'Difference'), clim=clim) p.add_time_browser(update_fields_function, t) p.show() @@ -264,18 +275,17 @@ def update_fields_function(i): # Account for the size of the source in the analytical formula fn_kdomain_finite_size = int_Fraunhofer_2D['gaussian'](R0_src) - ### -> Exact solution, At few points + # -> Exact solution, At few points x = points_out.T - r = np.linalg.norm(x - X0_src[np.newaxis,:], axis=1) - t = dt*np.arange(num_steps) + r = np.linalg.norm(x - X0_src[np.newaxis, :], axis=1) + t = dt * np.arange(num_steps) signals_exact = u_2D_SH_rt(r, src_t(t), rho.value, mu.value, dt, fn_kdomain_finite_size) # - fig, ax = plt.subplots(1,1) + fig, ax = plt.subplots(1, 1) ax.set_title('Signals at few points') for i in range(len(all_signals)): - ax.plot(t, all_signals[i,0,:], c='C'+str(i), ls='-') - ax.plot(t, signals_exact[i,:], c='C'+str(i), ls='--') + ax.plot(t, all_signals[i, 0, :], c=f'C{i}', ls='-') + ax.plot(t, signals_exact[i, :], c=f'C{i}', ls='--') ax.set_xlabel('Time') ax.legend(['FEM', 'analytical']) plt.show() - diff --git a/demo/weqnl_1D-PSV_Murnaghan_Pwave.py b/demo/weqnl_1D-PSV_Murnaghan_Pwave.py index 2d22599..655ef65 100644 --- a/demo/weqnl_1D-PSV_Murnaghan_Pwave.py +++ b/demo/weqnl_1D-PSV_Murnaghan_Pwave.py @@ -28,10 +28,10 @@ from mpi4py import MPI from petsc4py import PETSc -from elastodynamicsx.pde import material, BodyForce, boundarycondition, PDE, PDECONFIG +from elastodynamicsx.pde import material, boundarycondition, PDE, PDECONFIG from elastodynamicsx.solvers import TimeStepper from elastodynamicsx.plot import plotter -from elastodynamicsx.utils import spectral_element, spectral_quadrature, make_facet_tags, make_cell_tags, ParallelEvaluator +from elastodynamicsx.utils import spectral_element, spectral_quadrature, make_facet_tags, ParallelEvaluator # - # ### Set up a Spectral Element Method @@ -60,8 +60,8 @@ # define some tags tag_left, tag_right = 1, 2 all_tags = (tag_left, tag_right) -boundaries = [(tag_left , lambda x: np.isclose(x[0], 0 )),\ - (tag_right , lambda x: np.isclose(x[0], length))] +boundaries = [(tag_left , lambda x: np.isclose(x[0], 0 )), + (tag_right, lambda x: np.isclose(x[0], length))] facet_tags = make_facet_tags(domain, boundaries) # - @@ -96,8 +96,8 @@ # + T_N = fem.Constant(domain, default_scalar_type([0, 0])) # normal traction (Neumann boundary condition) Z_N, Z_T = mat.Z_N, mat.Z_T # P and S mechanical impedances -bc_l = boundarycondition((V, facet_tags, tag_left ), 'Neumann', T_N) -bc_rl = boundarycondition((V, facet_tags, (tag_left,tag_right) ), 'Dashpot', Z_N, Z_T) +bc_l = boundarycondition((V, facet_tags, tag_left), 'Neumann', T_N) +bc_rl = boundarycondition((V, facet_tags, (tag_left, tag_right)), 'Dashpot', Z_N, Z_T) bcs = [bc_l, bc_rl] # - @@ -111,17 +111,21 @@ T0 = 1 / f0 # period d0 = 5 * T0 # duration of source + def src_t(t): # source(t): Sine x Hann window - window = np.sin(np.pi*t / d0)**2 * (t < d0) * (t > 0) # Hann window - return np.sin(2*np.pi*f0 * t) * window + window = np.sin(np.pi * t / d0)**2 * (t < d0) * (t > 0) # Hann window + return np.sin(2 * np.pi * f0 * t) * window + # -> Space-Time function -p0 = default_scalar_type(7e-2) # max amplitude -F_0 = p0 * default_scalar_type([1, 1]) # source orientation +p_0 = default_scalar_type(7e-2) # max amplitude +F_0 = p_0 * default_scalar_type([1, 1]) # source orientation + def T_N_function(t): return src_t(t) * F_0 + if domain.comm.rank == 0: import matplotlib.pyplot as plt t = np.linspace(0, 1.1 * d0) @@ -187,17 +191,19 @@ def T_N_function(t): V.num_sub_spaces, num_steps)) # <- output stored here + # -> Define callbacks: will be called at the end of each iteration def cbck_storeAtPoints(i, out): if paraEval.nb_points_local > 0: - signals_local[:,:,i+1] = u_res.eval(paraEval.points_local, paraEval.cells_local) + signals_local[:, :, i+1] = u_res.eval(paraEval.points_local, paraEval.cells_local) + # enable live plotting enable_plot = True clim = 0.1 * np.amax(F_0) * np.array([0, 1]) -kwplot = { 'clim':clim, 'warp_factor':0.5 / np.amax(clim) } +kwplot = {'clim': clim, 'warp_factor': 0.5 / np.amax(clim)} if domain.comm.rank == 0 and enable_plot: - p = plotter(u_res, refresh_step=10, **kwplot) + p = plotter(u_res, refresh_step=10, window_size=[640, 480], **kwplot) if paraEval.nb_points_local > 0: # add points to live_plotter p.add_points(paraEval.points_local, render_points_as_spheres=True, point_size=12) @@ -215,6 +221,7 @@ def cbck_storeAtPoints(i, out): def cfst_updateSources(t): T_N.value = T_N_function(t) + # Run the big time loop! tStepper.solve(num_steps - 1, callfirsts=[cfst_updateSources], @@ -232,34 +239,39 @@ def cfst_updateSources(t): all_signals = paraEval.gather(signals_local, root=0) if domain.comm.rank == 0: - ### -> Exact (linear) solution, At few points + # -> Exact (linear) solution, At few points x = points_out.T - t = dt*np.arange(num_steps) - cL, cS = np.sqrt((lambda_.value+2*mu.value)/rho.value), np.sqrt(mu.value/rho.value) - resp_L = np.cumsum( src_t(t[np.newaxis,:]-x[:,0,np.newaxis]/cL) * F_0[0]/2/cL/rho.value, axis=1)*dt - resp_S = np.cumsum( src_t(t[np.newaxis,:]-x[:,0,np.newaxis]/cS) * F_0[1]/2/cS/rho.value, axis=1)*dt + t = dt * np.arange(num_steps) + cL = np.sqrt((lambda_.value + 2 * mu.value) / rho.value) + cS = np.sqrt(mu.value / rho.value) + resp_L = np.cumsum(src_t(t[np.newaxis, :] - x[:, 0, np.newaxis] / cL) * F_0[0] / 2 / cL / rho.value, axis=1) * dt + resp_S = np.cumsum(src_t(t[np.newaxis, :] - x[:, 0, np.newaxis] / cS) * F_0[1] / 2 / cS / rho.value, axis=1) * dt signals_linear = np.stack((resp_L, resp_S), axis=1) - # - f = np.fft.rfftfreq(t.size)/dt - fig, ax = plt.subplots(V.num_sub_spaces,2, sharex='col', sharey='none') - ax[0,0].set_title('Signals at few points') - ax[0,1].set_title('Spectra at few points') + + f = np.fft.rfftfreq(t.size) / dt + + fig, ax = plt.subplots(V.num_sub_spaces, 2, sharex='col', sharey='none') + ax[0, 0].set_title('Signals at few points') + ax[0, 1].set_title('Spectra at few points') for icomp, cax in enumerate(ax): for i in range(len(all_signals)): - cax[0].text(0.02,0.97, 'U'+['x','y','z'][icomp], ha='left', va='top', transform=cax[0].transAxes) - cax[0].plot(t, all_signals[i,icomp,:], c='C'+str(i), ls='-' ) # FEM - cax[0].plot(t, signals_linear[i,icomp,:], c='C'+str(i), ls='--') # exact linear + cax[0].text(0.02, 0.97, 'U' + ['x', 'y', 'z'][icomp], ha='left', va='top', transform=cax[0].transAxes) + cax[0].plot(t, all_signals[i, icomp, :], c=f'C{i}', ls='-') # FEM + cax[0].plot(t, signals_linear[i, icomp, :], c=f'C{i}', ls='--') # exact linear # - cax[1].plot(f, np.abs(np.fft.rfft(all_signals[i,icomp,:])), c='C'+str(i), ls='-' ) # FEM - cax[1].plot(f, np.abs(np.fft.rfft(signals_linear[i,icomp,:])), c='C'+str(i), ls='--') # exact linear - specX = np.abs(np.fft.rfft(all_signals[i,0,:])) - specX_f, specX_2f = specX[np.argmin(np.abs(f-1/T0))], specX[np.argmin(np.abs(f-2/T0))] - ax[0, 1].annotate('2nd harmonic', xy=(2/T0, 1.2*specX_2f), xytext=(2.1/T0, 0.4*specX_f), arrowprops=dict(facecolor='black', shrink=0.05)) - ax[-1,0].set_xlabel('Time') - ax[-1,1].set_xlabel('Frequency') - ax[-1,1].set_xlim(-0.1*1/T0, 3.8*1/T0) - ax[-1,1].set_xticks(np.arange(4)/T0, ['0', 'f', '2f', '3f']) + cax[1].plot(f, np.abs(np.fft.rfft(all_signals[i, icomp, :])), c=f'C{i}', ls='-') # FEM + cax[1].plot(f, np.abs(np.fft.rfft(signals_linear[i, icomp, :])), c=f'C{i}', ls='--') # exact linear + + specX = np.abs(np.fft.rfft(all_signals[i, 0, :])) + specX_f = specX[np.argmin(np.abs(f - 1 / T0))] + specX_2f = specX[np.argmin(np.abs(f - 2 / T0))] + ax[0, 1].annotate('2nd harmonic', xy=(2 / T0, 1.2 * specX_2f), + xytext=(2.1 / T0, 0.4 * specX_f), + arrowprops=dict(facecolor='black', shrink=0.05)) + ax[-1, 0].set_xlabel('Time') + ax[-1, 1].set_xlabel('Frequency') + ax[-1, 1].set_xlim(-0.1 / T0, 3.8 / T0) + ax[-1, 1].set_xticks(np.arange(4) / T0, ['0', 'f', '2f', '3f']) plt.show() # # ----------------------------------------------------- - From 74b261db77b0c37ae901952bb00cf9eb7c2caf82 Mon Sep 17 00:00:00 2001 From: Pierric Mora Date: Tue, 23 Apr 2024 12:05:13 +0200 Subject: [PATCH 12/12] v0.2.3 --- elastodynamicsx/_version.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/elastodynamicsx/_version.py b/elastodynamicsx/_version.py index dd197b1..f4a41f3 100644 --- a/elastodynamicsx/_version.py +++ b/elastodynamicsx/_version.py @@ -4,4 +4,4 @@ # # SPDX-License-Identifier: MIT -__version__ = '0.2.2' +__version__ = '0.2.3'