Wave equation simulation with dolfinx

so I came across this code for the simulation of acoustic wave equation.
fenics-scripts/awefem/awefem.py at master · cako/fenics-scripts · GitHub

has anyone rewritten this for dolfinx? I’m to new to fenics and just needed the simulation. I tried to install legacy versions with docker but it’s pulling nothing.

You can get legacy fenics docker images at: Package fenics-gmsh · GitHub

Looking at the code, it shouldn’t be too hard to write for DOLFINx, as you for instance can get point sources through scifem: Point sources in DOLFINx — scifem

1 Like

thanks for the quick reply. I couldn’t install scifem it had some issues with version of CMAKE. but here I tried to implement, can you check it out as I see no output:

from vedo import download
from vedo.dolfin import plot

import importlib.util
from mpi4py import MPI
import numpy as np
import pyvista as pv
from ufl import *
from dolfinx import *
from dolfinx.io import XDMFFile
from dolfinx.plot import vtk_mesh
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, set_bc, LinearProblem
import basix.ufl
if importlib.util.find_spec("petsc4py") is not None:
    import dolfinx

    if not dolfinx.has_petsc:
        print("This demo requires DOLFINx to be compiled with PETSc enabled.")
        exit(0)
    from petsc4py import PETSc
    from petsc4py.PETSc import ScalarType,KSP,PC,InsertMode  # type: ignore
else:
    print("This demo requires petsc4py.")
    exit(0)

def ricker_source(t, f=40):
    t -= 2 / f
    return (1 - 2 * (np.pi*f*t)**2) * np.exp(-(np.pi*f*t)**2)

def sine_source(t, f=40):
    return np.sin(2 * np.pi*f*t)


def awefem(msh, t, source_loc=None):
    # Function space
    V = fem.functionspace(msh, ("Lagrange", 1))
    u_bc_value = fem.Constant(msh, ScalarType(0))
    boundary_facets = mesh.locate_entities_boundary(
        msh, 
        msh.topology.dim - 1,
        lambda x: np.full(x.shape[1], True)
        )
    dofs = fem.locate_dofs_topological(V, msh.topology.dim - 1, boundary_facets)
    bc = fem.dirichletbc(u_bc_value, dofs, V)
    

    # Trial and test functions
    u = TrialFunction(V)
    v = TestFunction(V)

    # Discretization
    c = 6
    dt = t[1] - t[0]
    u0 = fem.Function(V)  # u0 = uN-1
    u1 = fem.Function(V)  # u1 = uN1

    # Variational formulation
    # F = (u - 2 * u1 + u0) * v * dx + (dt * c) ** 2 * dot(grad(u + 2 * u1 + u0) / 4, grad(v) ) * dx
    F = (u - 2 * u1 + u0) * conj(v) * dx + (dt * c) ** 2 * dot(grad(u + 2 * u1 + u0) / 4, grad(conj(v))) * dx
    a, L = fem.form(lhs(F)), fem.form(rhs(F))

    # Solver
    A = assemble_matrix(a)
    A.assemble()
    b = fem.assemble_vector(L)
    fem.apply_lifting(b.array, [a], bcs=[[bc]])
    bc.set(b.array)
    solver = PETSc.KSP().create(msh.comm)
    solver.setType(PETSc.KSP.Type.PREONLY)  
    solver.getPC().setType(PETSc.PC.Type.LU) 
    solver.setOperators(A)

    # Solution
    u = fem.Function(V)  # uN+1

    # Source
    if source_loc is None:
        mesh_center = np.mean(msh.geometry.x, axis=0)
        source_loc = tuple(mesh_center)
    else:
        source_loc = tuple(source_loc)
    source_loc = np.array([source_loc[0], source_loc[1], 0.0])
    # with XDMFFile(msh.comm, "square/data/square.xdmf", "w") as xdmf:
    #     xdmf.write_mesh(msh)
    
    # Set up the PyVista plotter for animation
    topology, cells, geometry = vtk_mesh(msh)
    grid = pv.UnstructuredGrid(topology, cells, geometry)
    u_values = u.x.array  # Initial data for the solution
    grid["u"] = u_values

    plotter = pv.Plotter()
    plotter.add_mesh(grid, scalars="u", clim=[-1, 1], cmap="seismic")
    plotter.show_axes()
    plotter.add_scalar_bar("u", title="Amplitude")
    plotter.show(interactive_update=True, auto_close=False)
    
    # Time stepping
    for i, t_ in enumerate(t[1:]):
        b = assemble_vector(L)
        source_dof = fem.locate_dofs_geometrical(
            V, lambda x: np.isclose(x.T, source_loc).all(axis=1)
            )
        
        if len(source_dof) == 0:
            raise ValueError(f"No DOFs found near source location {source_loc}")

        point_source_value = ricker_source(t_) * dt**2
        b.setValues(source_dof, [point_source_value], addv=PETSc.InsertMode.ADD)
        b.assemble()

        solver.solve(b, u.x.petsc_vec)
        u0.x.array[:] = u1.x.array
        u1.x.array[:] = u.x.array
        
        # xdmf.write_function(u, t_)
        u_values[:] = u.x.array
        grid["u"] = u_values
        plotter.update_scalars(u_values, render=True)
if __name__ == "__main__":

    ot, dt, nt = 0, 5e-4, 1001
    t = ot + np.arange(nt) * dt

    # WIP: couldn't download
    # print("Computing wavefields over dolfin mesh")
    # fpath = download("https://vedo.embl.es/examples/data/dolfin_fine.xml") 
    # msh = Mesh(fpath)
    
    # N0=7
    # N = N0 * 2**5
    # comm = MPI.COMM_WORLD
    # msh = mesh.create_unit_square(comm, N, N, mesh.CellType.triangle)

    msh = mesh.create_rectangle(
        comm=MPI.COMM_WORLD,
        points=((0.0, 0.0), (2.0, 1.0)),
        n=(100, 100),
        cell_type=mesh.CellType.triangle,
        )
    awefem(msh, t, source_loc=(1, 1))

Could you be more specific? What was the error message.

Note that there is a point source implementation at:
Point Sources (Redux) - #9 by francesco-ballarin Which is what I based the scifem implementation of

sure here’s the output:

error: subprocess-exited-with-error

  × Building wheel for scifem (pyproject.toml) did not run successfully.
  │ exit code: 1
  ╰─> [35 lines of output]
      *** scikit-build-core 0.10.7 using CMake 3.24.4 (wheel)
      *** Configuring CMake...
      2024-11-19 13:01:50,383 - scikit_build_core - WARNING - Unsupported CMAKE_ARGS ignored: -DCMAKE_BUILD_TYPE=Release
      2024-11-19 13:01:50,457 - scikit_build_core - WARNING - Unsupported CMAKE_ARGS ignored: -DCMAKE_BUILD_TYPE=Release
      loading initial cache file /var/folders/k9/zmdp8k115vqdx4bgg6dn7dlh0000gn/T/tmpcv0h7r_d/build/CMakeInit.txt
      -- The C compiler identification is Clang 18.1.8
      -- The CXX compiler identification is Clang 18.1.8
      -- Detecting C compiler ABI info
      -- Detecting C compiler ABI info - done
      -- Check for working C compiler: /Users/mjazbay/miniconda3/envs/fenicsx-env/bin/arm64-apple-darwin20.0.0-clang - skipped
      -- Detecting C compile features
      -- Detecting C compile features - done
      -- Detecting CXX compiler ABI info
      -- Detecting CXX compiler ABI info - done
      -- Check for working CXX compiler: /Users/mjazbay/miniconda3/envs/fenicsx-env/bin/clang++ - skipped
      -- Detecting CXX compile features
      -- Detecting CXX compile features - done
      -- Found Basix at /Users/mjazbay/miniconda3/envs/fenicsx-env/lib/cmake/basix
      -- Could NOT find MPI_C (missing: MPI_C_WORKS)
      -- Could NOT find MPI_CXX (missing: MPI_CXX_WORKS)
      CMake Error at /opt/local/share/cmake-3.24/Modules/FindPackageHandleStandardArgs.cmake:230 (message):
        Could NOT find MPI (missing: MPI_C_FOUND MPI_CXX_FOUND)
      Call Stack (most recent call first):
        /opt/local/share/cmake-3.24/Modules/FindPackageHandleStandardArgs.cmake:594 (_FPHSA_FAILURE_MESSAGE)
        /opt/local/share/cmake-3.24/Modules/FindMPI.cmake:1835 (find_package_handle_standard_args)
        /opt/local/share/cmake-3.24/Modules/CMakeFindDependencyMacro.cmake:47 (find_package)
        /Users/mjazbay/miniconda3/envs/fenicsx-env/lib/cmake/dolfinx/DOLFINXConfig.cmake:35 (find_dependency)
        CMakeLists.txt:35 (find_package)


      -- Configuring incomplete, errors occurred!
      See also "/var/folders/k9/zmdp8k115vqdx4bgg6dn7dlh0000gn/T/tmpcv0h7r_d/build/CMakeFiles/CMakeOutput.log".
      See also "/var/folders/k9/zmdp8k115vqdx4bgg6dn7dlh0000gn/T/tmpcv0h7r_d/build/CMakeFiles/CMakeError.log".

      *** CMake configuration failed
      [end of output]

  note: This error originates from a subprocess, and is likely not a problem with pip.
  ERROR: Failed building wheel for scifem
ERROR: ERROR: Failed to build installable wheels for some pyproject.toml based projects (scifem)

Did you use conda install -c conda-forge scifem?

no, I used python3 -m pip install scifem --no-build-isolation and also tried python3 -m pip install --no-build-isolation git+https://github.com/scientificcomputing/scifem.git

Since you are using conda I would use the conda installation of scifem, to ensure that everything is compatible: Scifem | Anaconda.org

It’s also unccessfull using conda

Channels:
 - defaults
 - conda-forge
Platform: osx-arm64
Collecting package metadata (repodata.json): done
Solving environment: | warning  libmamba Added empty dependency for problem type SOLVER_RULE_UPDATE
failed

LibMambaUnsatisfiableError: Encountered problems while solving:
  - nothing provides _python_rc needed by python-3.12.0rc3-rc3_h47c9636_1_cpython

Could not solve for environment specs
The following packages are incompatible
├─ pin-1 is installable and it requires
│  └─ python 3.13.* , which can be installed;
└─ scifem is not installable because there are no viable options
   ├─ scifem [0.2.5|0.2.7|0.2.8] would require
   │  └─ python 3.10.* *_cpython, which conflicts with any installable versions previously reported;
   ├─ scifem [0.2.5|0.2.7|0.2.8] would require
   │  └─ python 3.11.* *_cpython, which conflicts with any installable versions previously reported;
   └─ scifem [0.2.5|0.2.7|0.2.8] would require
      └─ python 3.12.* *_cpython but there are no viable options
         ├─ python [3.12.0|3.12.1|...|3.12.7] conflicts with any installable versions previously reported;
         └─ python 3.12.0rc3 would require
            └─ _python_rc, which does not exist (perhaps a missing channel).

Please try with a new environment using python3.12. We haven’t had time to build Python 3.13 binaries yet.

okay, thanks that worked for me with python 3.12.2. I still couldn’t run that example with point_source.py. My fenics-dolfinx 0.9.0

Traceback (most recent call last):
  File "/Users/mjazbay/fenicsx/awefem/point_source.py", line 98, in <module>
    dolfinx.fem.petsc.apply_lifting(b.vector, [a_compiled], [[bc]])
                                    ^^^^^^^^
AttributeError: 'Function' object has no attribute 'vector'

I have not seen the specific model you are following but b.vector is deprecated as much as I remember. Try b.x.petsc_vec instead

thanks, this fixed the above problem.

is there images of legacy versions before the latest? I’m trying to run this legacy script but it crashes for the latest version.

There are plenty of releases:

yes there’s plenty, but all of them:

fenics-dijitso    2019.2.0.dev0
fenics-dolfin     2019.2.0.dev0
fenics-ffc        2019.2.0.dev0
fenics-fiat       2019.2.0.dev0
fenics-ufl-legacy 2022.3.0

I guess your error is related to ufl/ufl_legacy?

There was a switch around April 2023, so an image prior to that should resolve it.

We do not have the bandwidth to build images of releases that are more than 5 years old.