How to get the solution value at the boundary or any give points?

General
1

Hello everyone. I am solving a set of coupled 1d-PDEs,

\begin{cases} \frac{\partial c_1}{\partial t} &= \frac{\partial}{\partial x}f_1(c_1,c_2,...,\partial c_1/\partial x, \partial c_2/\partial x,...) \\ \frac{\partial c_2}{\partial t} &= \frac{\partial}{\partial x}f_2(c_1,c_2,...,\partial c_1/\partial x, \partial c_2/\partial x,...) \end{cases}

the left boundary is the flux \vec{f}(c_1, c_2, ...), it’s a numerical function, so I can’t write it in the variational form. How can I get the \{c_i\} and \{\partial c_i/\partial x\} at the left boundary? A more general question is, how can I get them at arbitrary points on the domain?

I find a topic in Calculate the value of a solution on mesh points or in any point - General - FEniCS Project
But it does not works. The error occurs at

bb_tree = geometry.BoundingBoxTree(mesh, mesh.topology.dim)

python3 test_eval.py 
Traceback (most recent call last):
  File "/mnt/g/workdir/MKM/NO_mkm/test_eval.py", line 46, in <module>
    bb_tree = geometry.BoundingBoxTree(mesh, mesh.topology.dim)
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
TypeError: BoundingBoxTree.__init__() takes 2 positional arguments but 3 were given

And I don’t understand these codes at all because the lack of comments. For example, the functions such as

exterior_facet_indices
tabulate_dof_coordinates
compute_collisions
compute_colliding_cells
...
2

There is a test case, a vector-valued problem on (0, 1)

\begin{cases} \frac{\partial u_1}{\partial t} &= \frac{\partial^2 u_1}{\partial x^2} - \exp(u_1 - u_2) + \exp(-(u_1 - u_2)) \\ \frac{\partial u_2}{\partial t} &= \frac{\partial^2 u_2}{\partial x^2} + \exp(u_1 - u_2) - \exp(-(u_1 - u_2)) \\ \end{cases}

the boundary condition is

\begin{cases} \frac{\partial u_1}{\partial x}(0, t) &= \sin(u_1(x = 0)) \\ \frac{\partial u_2}{\partial x}(0, t) &= -0.1u_1(x=0)u_2(x=0) \end{cases} \quad \begin{cases} u_1(1, t) = 1 \\ u_2(1, t) = 0 \end{cases} 
# -*- coding:utf-8 -*-
import ufl 
import dolfinx
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from petsc4py import PETSc
from mpi4py import MPI
import numpy as np 

t_ini = 0
t_fin = 1.1
n_step = 11
dt = (t_fin - t_ini) / n_step

nx = 64
L = 1.0
rank = MPI.COMM_WORLD.Get_rank()

domain = dolfinx.mesh.create_interval(MPI.COMM_WORLD, nx, points=(0.0, L))
V = dolfinx.fem.FunctionSpace(domain, ufl.VectorElement("Lagrange", domain.ufl_cell(), 1, dim=2))

fdim = domain.topology.dim - 1

boundary_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(L, x[0]))

bc1 = dolfinx.fem.dirichletbc(PETSc.ScalarType(1.0), 
                              dolfinx.fem.locate_dofs_topological(V.sub(0), fdim, boundary_facets), 
                              V.sub(0))
bc2 = dolfinx.fem.dirichletbc(PETSc.ScalarType(0.0), 
                              dolfinx.fem.locate_dofs_topological(V.sub(1), fdim, boundary_facets), 
                              V.sub(1))


u_n = dolfinx.fem.Function(V)
u = dolfinx.fem.Function(V)
u_n1, u_n2 = ufl.split(u_n)
u1, u2 = ufl.split(u)

u_n.sub(0).interpolate(lambda x: np.ones(x[0].shape))
u_n.sub(1).interpolate(lambda x: np.zeros(x[0].shape))
u_n.x.scatter_forward()

v1, v2 = ufl.TestFunction(V)
G = -ufl.exp(u1-u2)+ufl.exp(-(u1-u2)) 
F1 = u1*v1*ufl.dx + dt*ufl.dot(ufl.grad(u1), ufl.grad(v1))*ufl.dx + dt*ufl.sin(u1)*v1*ufl.ds - dt*G*v1*ufl.dx - u_n1*v1*ufl.dx
F2 = u2*v2*ufl.dx + dt*ufl.dot(ufl.grad(u2), ufl.grad(v2))*ufl.dx + dt*(-0.1*u1*u2)*v2*ufl.ds + dt*G*v2*ufl.dx - u_n2*v2*ufl.dx
F = F1 + F2


problem = NonlinearProblem(F, u, bcs=[bc1, bc2])
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "residual"
solver.rtol = 1e-6

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

file = dolfinx.io.XDMFFile(MPI.COMM_WORLD, "diffu_test3/u_tot.xdmf", "w")
file.write_mesh(domain)

t = t_ini 
while t < t_fin:
    n, converged = solver.solve(u)
    if (rank == 0):
        print(f"Number of interations: {n:d}")
    u_n.x.array[:] = u.x.array
    file.write_function(u, t)
    t += dt
    
file.close()
3

See
https://jsdokken.com/dolfinx-tutorial/chapter1/membrane_code.html#making-curve-plots-throughout-the-domain
for updated syntax regarding boundingboxes

1 Like
4

Here is a working example which evaluate the solution at the left boundary (0.0, 0.0, 0.0) and broad cast the value to all the process.

cells = []
points_on_proc = []
bb_tree = dolfinx.geometry.bb_tree(domain, domain.topology.dim)
points = np.array([[0.0, 0.0, 0.0]])
cell_candidates = dolfinx.geometry.compute_collisions_points(bb_tree, points)
colliding_cells = dolfinx.geometry.compute_colliding_cells(domain, cell_candidates, points)
for i, point in enumerate(points):
    if len(colliding_cells.links(i)) > 0:
        points_on_proc.append(point)
        cells.append(colliding_cells.links(i)[0])
points_on_proc = np.array(points_on_proc, dtype=np.float64)

non_empty_process = MPI.COMM_WORLD.gather(len(points_on_proc), root=0)
if rank == 0:
    root_rank = np.argmax(non_empty_process)
else:
    root_rank = None
root_rank = MPI.COMM_WORLD.bcast(root_rank, root=0)

t = t_ini 
while t < t_fin:
    n, converged = solver.solve(u)
    
    u_values = u.eval(points_on_proc, cells)
    u_values = MPI.COMM_WORLD.bcast(u_values, root=root_rank)
    print("rank = ", rank, "u_values = ", u_values)

    u_n.x.array[:] = u.x.array
    file.write_function(u, t)
    t += dt
5

Thanks for your answer, my problem is solved.