How to implement point source in elastic wave propagation

[tengfeideng]

Hello, guys,
I want to ask how to implement point source in Firedrake. And I found a reference code about ithttps://www.firedrakeproject.org/demos/higher_order_mass_lumping.py.html
def RickerWavelet(t, freq, amp=1.0):
t = t - (np.sqrt(6.0) / (np.pi * freq))
return amp * (1.0 - (1.0 / 2.0) * (2.0 * np.pi * freq) * (2.0 * np.pi * freq) * t * t)

def delta_expr(x0, x, y, sigma_x=2000.0):
sigma_x = Constant(sigma_x)
return exp(-sigma_x * ((x - x0[0]) ** 2 + (y - x0[1]) ** 2))
But It is defined in a FunctionSpace. Does anybody know how to define a PointSource in VectorFunctionSpace? I want to simulate an explosive source.

Thanks a lot.

Tengfei

[lrtfm]

The return value of delta_expr is a scalar ufl expression. Do you want to get a vector expression?
If that is what you want, you can do this:

def delta_expr_2d(x0, x, y, sigma_x=2000.0):
  sigma_x = Constant(sigma_x)
  exp1 = exp(-sigma_x * ((x - x0[0]) ** 2 + (y - x0[1]) ** 2))  
  exp2 = exp(-sigma_x * ((x - x0[0]) ** 2 + (y - x0[1]) ** 2))  
  return as_vector( (exp1, exp2) )

[PG20R]

Hi. I tried to implement a simple scalar delta function as in your discussion. I applied my delta to three test cases as shown below and it does not work at all. My test code:

from firedrake import *

test_mesh = UnitDiskMesh(5)
print('num_vertices', test_mesh.num_vertices())

def delta_expr(x0, x,y, sigma_x=2000.0):
    sigma_x = Constant(sigma_x)
    return exp(-sigma_x * ((x-x0[0])**2+(y-x0[1])**2))
    
x,y = SpatialCoordinate(test_mesh)

source = Constant([0.0,0.0])
delta = delta_expr(source,x,y)
test1 = (x*y)**2 - Constant(1.0)
test2 = x+y
test3 = cos(x*y)
print('Test1', assemble(delta*test1*dx))  #expect -1.0
print('Test2', assemble(delta*test2*dx))   #expect 0.0
print('Test3', assemble(delta*test3*dx))   #expect 1.0

What is going wrong here? Do I have to adjust sigma_x to my mesh size? If yes, how to do this? Or how else could I define a simple point source in Firedrake?

[lrtfm]

Here, as you use an approximation to get the Dirac delta function, you need a proper coefficient for the return expression, see https://en.wikipedia.org/wiki/Dirac_delta_function for some approximations. However, those approximations may only work when your point is in the domain, not on the boundary. On the boundary, the coefficients depend on the angle of the point.

Here, another implementation is using quadrature rule:

from firedrake import *

from firedrake.petsc import PETSc
from pyop2 import op2
from pyop2.datatypes import ScalarType
from mpi4py import MPI
import finat
import numpy as np

set_level(CRITICAL) # Disbale warnings


def delta_expr(m, x0):
    """Make dirac function at point

    Args:
        m: mesh
        x0: source point
    Return:
        delta function

    Example:
        delta = delta_expr(m, x0)
        f = Function(V)
        f_x0 = assemble(delta(f))
    """
    V = FunctionSpace(m, 'DG', 0)
    cell_marker = Function(V, name='cell_marker', dtype=ScalarType)
    qrule = finat.quadrature.make_quadrature(V.finat_element.cell, 0)
    cell, X = m.locate_cell_and_reference_coordinate(x0, tolerance=1e-6)

    # c = 0 if X is None else 1
    n_cell_local = len(cell_marker.dat.data)
    if X is not None and cell < n_cell_local:
        c = 1
    else:
        c = 0

    comm = m.comm
    s = comm.size - comm.rank
    n = comm.allreduce(int(s*c), op=MPI.MAX)

    if n == 0:
        raise BaseException("Points not found!")

    k = int(comm.size - n) # get the lower rank which include the point x0

    if c == 1 and comm.rank == k:
        X[X<0] = 0
        X[X>1] = 1
        cell_marker.dat.data[cell] = 1
        comm.bcast(X, root=k)
    else:
        cell_marker.dat.data[:] = 0 # we must set this otherwise the process will hangup
        X = comm.bcast(None, root=k)


    cell_marker.dat.global_to_local_begin(op2.READ)
    cell_marker.dat.global_to_local_end(op2.READ)

    qrule.point_set.points[0] = X
    qrule.weights[0] = qrule.weights[0]/np.real(assemble(cell_marker*dx))

    return lambda f: f*cell_marker*dx(rule=qrule)


if __name__ == '__main__':
    test_mesh = UnitDiskMesh(5)
    v_num = test_mesh.topology_dm.getVertexNumbering()
    m_local = max(v_num.array)
    num_vertices = test_mesh.comm.allreduce(m_local, MPI.MAX)

    PETSc.Sys.Print('num_vertices', num_vertices)

    x = SpatialCoordinate(test_mesh)

    source = Constant([0.0,0.0])
    delta = delta_expr(test_mesh, source)
    test1 = (x[0]*x[1])**2 - Constant(1.0)
    test2 = x[0]+x[1]
    test3 = cos(x[0]*x[1])
    PETSc.Sys.Print('Test1', assemble(delta(test1)))  #expect -1.0
    PETSc.Sys.Print('Test2', assemble(delta(test2)))   #expect 0.0
    PETSc.Sys.Print('Test3', assemble(delta(test3)))   #expect 1.0

$ mpiexec -n 4 python dirac.py
num_vertices 4224
Test1 -0.9999999999999989
Test2 -4.661211986717451e-17
Test3 0.9999999999999989
$ mpiexec -n 1 python dirac.py
num_vertices 4224
Test1 -0.9999999999999989
Test2 7.703719777548935e-34
Test3 0.9999999999999989

[colinjcotter]

What do you get? You didn't show it but I assume you are getting different numbers. It is expected that you won't get the exact answer, because Firedrake doesn't do exact integration, but numerical quadrature (of course). And you have a 5x5 mesh so you have very few triangles. The quadrature rules are such that the right answer will be achieved in the limit of large numbers of triangles.

…

________________________________________ From: PG20R ***@***.***> Sent: 22 August 2022 09:25 To: firedrakeproject/firedrake Cc: Subscribed Subject: Re: [firedrakeproject/firedrake] How to implement point source in elastic wave propagation (Discussion #2510) Hi. I tried to implement a simple scalar delta function as in your discussion. I applied my delta to three test cases as shown below and it does not work at all. My test code: from firedrake import * test_mesh = UnitDiskMesh(5) print('num_vertices', test_mesh.num_vertices()) def delta_expr(x0, x,y, sigma_x=2000.0): sigma_x = Constant(sigma_x) return exp(-sigma_x * ((x-x0[0])**2+(y-x0[1])**2)) x,y = SpatialCoordinate(test_mesh) source = Constant([0.0,0.0]) delta = delta_expr(source,x,y) test1 = (x*y)**2 - Constant(1.0) test2 = x+y test3 = cos(x*y) print('Test1', assemble(delta*test1*dx)) #expect -1.0 print('Test2', assemble(delta*test2*dx)) #expect 0.0 print('Test3', assemble(delta*test3*dx)) #expect 1.0 What is going wrong here? Do I gave to adjust sigma_x to my mesh size? If yes, how to this? Or how else could I define a simple point source in Firedrake? — Reply to this email directly, view it on GitHub<#2510 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABOSV4UOWKQJHOYSEI2YQD3V2M2ONANCNFSM55EL5SLQ>. You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>