Green's function for advetion diffusion operator in 2D

General
1

Hi everyone,

I’m looking for resources on computing the Green’s function for a 2D advection-diffusion equation using FEniCSx. Ideally, I’d like an example that demonstrates this process.

Additionally, I’m interested in knowing how to implement a delta function as a forcing term in FEniCSx. While I understand this might be a general question, any pointers on achieving this or relevant resources would be greatly appreciated.

Thanks in advance!

2

For part two of your question see

1 Like
3

Thank you dr. @dokken for your prompt response. I will study the code thoroughly and will follow up if I have any issues. I hope I am not bothering.

4

Hii dr. @dokken I tried to run the code for point source evaluation as it is: it is giving following error, I wonder what is going wrong:

# Author: Jørgen S. Dokken
#
# Create a point source

import dolfinx
from mpi4py import MPI
import numpy as np
import ufl


def compute_cell_contributions(V, points):
    # Determine what process owns a point and what cells it lies within
    _, _, owning_points, cells = dolfinx.cpp.geometry.determine_point_ownership(
        mesh._cpp_object, points, 1e-6)
    owning_points = np.asarray(owning_points).reshape(-1, 3)

    # Pull owning points back to reference cell
    mesh_nodes = mesh.geometry.x
    cmap = mesh.geometry.cmaps[0]
    ref_x = np.zeros((len(cells), mesh.geometry.dim),
                     dtype=mesh.geometry.x.dtype)
    for i, (point, cell) in enumerate(zip(owning_points, cells)):
        geom_dofs = mesh.geometry.dofmap[cell]
        ref_x[i] = cmap.pull_back(point.reshape(-1, 3), mesh_nodes[geom_dofs])

    # Create expression evaluating a trial function (i.e. just the basis function)
    u = ufl.TrialFunction(V)
    num_dofs = V.dofmap.dof_layout.num_dofs * V.dofmap.bs
    if len(cells) > 0:
        # NOTE: Expression lives on only this communicator rank
        expr = dolfinx.fem.Expression(u, ref_x, comm=MPI.COMM_SELF)
        values = expr.eval(mesh, np.asarray(cells))

        # Strip out basis function values per cell
        basis_values = values[:num_dofs:num_dofs*len(cells)]
    else:
        basis_values = np.zeros(
            (0, num_dofs), dtype=dolfinx.default_scalar_type)
    return cells, basis_values


mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

if mesh.comm.rank == 0:
    points = np.array([[0.1, 0.2, 0],
                      [0.91, 0.93, 0]], dtype=mesh.geometry.x.dtype)
else:
    points = np.zeros((0, 3), dtype=mesh.geometry.x.dtype)

V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))

cells, basis_values = compute_cell_contributions(V, points)
print(cells, basis_values)

Error

AttributeError                            Traceback (most recent call last)
<ipython-input-2-1a2a7fd3c87b> in <cell line: 52>()
     50 V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
     51 
---> 52 cells, basis_values = compute_cell_contributions(V, points)
     53 print(cells, basis_values)

<ipython-input-2-1a2a7fd3c87b> in compute_cell_contributions(V, points)
     17     # Pull owning points back to reference cell
     18     mesh_nodes = mesh.geometry.x
---> 19     cmap = mesh.geometry.cmaps[0]
     20     ref_x = np.zeros((len(cells), mesh.geometry.dim),
     21                      dtype=mesh.geometry.x.dtype)

AttributeError: 'Geometry_float64' object has no attribute 'cmaps'
5

change the offending line to
cmap = mesh.geometry.cmap

1 Like
6

Thank you so much dr. @dokken.