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