# PointSource in diffusion reaction problem

Hi there,

I am trying to model reaction-diffusion of biological cell signals, with the cells acting as point sources as they secrete the compound.

The code below runs but my approximation at point sources (supply term) is very clumsy. I have looked into posts such as PointSource in dolfinx - dolfinx - FEniCS Project and How to construct a sparse matrix used for calculating function values at multiple points - dolfinx - FEniCS Project but I don’t see how to integrate the true point source definition into my code.

Would you be able to point me in the right direction? The code has to work in parallel.

Also, would I be able to use the same approach to model point sinks ? (cell consumption)

Many thanks!

MWE:

``````L = 1
W = 0.2
diffusion_time_step=0.01
end_time=1
simu_start_time=0
max_dif_tstep=0.2

import numpy as np
import ufl

from mpi4py import MPI
from petsc4py.PETSc import ScalarType

from dolfinx import mesh, fem, plot, io, nls
from dolfinx.fem import FunctionSpace,VectorFunctionSpace, Function,TensorFunctionSpace, Constant
from ufl import TestFunction, TrialFunction, dx, inner,grad, Measure,derivative,system
from petsc4py import PETSc
import time

startT = time.time()

domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0,0,0]), np.array([L, W, W])],
[30,30,30], cell_type=mesh.CellType.hexahedron)

V=VectorFunctionSpace(domain, ("CG", 1), dim=3)

p = TestFunction(V)
prc_prev = Function(V)
prc = Function(V)

dt=Constant(domain, diffusion_time_step)

def set_time_step(dt,converged,n):
if converged and n<5:
new_dt=min(max_dif_tstep,2*dt)
else: new_dt=dt
return new_dt

D_p = Constant(domain,2.4e3) # Diffusion coefficient for PTHrP
D_pr = Constant(domain,2.4e-15) # Diffusion coefficient for PTHrP
kap = Constant(domain,1e5) # Binding association rate for PTHrP-receptor
kdp = Constant(domain,1.86e-3)  # Disassociation rate for PTHrP-receptor
psi_p = Constant(domain,0.01) # Decay rate for PTHrP
psi_pr = Constant(domain,1e-50) # Decay rate for PTHrP receptor
psi_cp = Constant(domain,1e-15) # Decay rate for PTHrP-receptor Complex

xdmf = io.XDMFFile(domain.comm, "MiniExampleCustN.xdmf", "w")
xdmf.write_mesh(domain)

# Supply term
def sup_func(x):
return np.exp(-100*(x[0]%0.1+x[1]%0.1+x[2]%0.1))
p_i = Function(V.sub(0).collapse()[0])
p_i.interpolate(sup_func)
p_sup = p_i

# PTHrP receptor presentation
def rec_func(x):
return np.exp(-10*(x[0]%0.15+x[1]%0.15+x[2]%0.15))
rp_i = Function(V.sub(0).collapse()[0])
rp_i.interpolate(rec_func)

subspace,submap=V.sub(1).collapse()
prc.x.array[submap]=rp_i.x.array
prc_prev.x.array[submap]=rp_i.x.array

# Variational problem
F = (
(- (prc[0] - prc_prev[0])*p[0]
- (prc[1] - prc_prev[1])*p[1]
- (prc[2] - prc_prev[2])*p[2])*dx()
+ (+ dt*p_sup*p[0])*dx()
+ (+ dt*(- kap*prc_prev[0]*prc[1]*p[0]
- kap*prc[0]*prc_prev[1]*p[1]
+ kap*prc_prev[0]*prc_prev[1]*p[2]
+ kdp*prc_prev[2]*p[0]
+ kdp*prc_prev[2]*p[1]
- kdp*prc[2]*p[2]))*dx()
+ (+ dt*(- psi_p*prc[0]*p[0]
- psi_pr*prc[1]*p[1]
- psi_cp*prc[2]*p[2]))*dx()
)

# Create boundary condition
tdim = V.mesh.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(V.mesh.topology)
fixed_dofs_0 = fem.locate_dofs_topological(V.sub(0), fdim, boundary_facets)
fixed_dofs_1 = fem.locate_dofs_topological(V.sub(1), fdim, boundary_facets)
fixed_dofs_2 = fem.locate_dofs_topological(V.sub(2), fdim, boundary_facets)
bc0 = fem.dirichletbc(PETSc.ScalarType(0), fixed_dofs_0, V.sub(0))
bc1 = fem.dirichletbc(PETSc.ScalarType(0), fixed_dofs_1, V.sub(1))
bc2 = fem.dirichletbc(PETSc.ScalarType(0), fixed_dofs_2, V.sub(2))
bc=[bc0,bc1,bc2]

# Define problem and custom Newton solver

solver = PETSc.KSP().create(domain.comm)
opts = PETSc.Options()
option_prefix =solver.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "hypre"
# opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
solver.setFromOptions()

residual = fem.form(F)
J = ufl.derivative(F, prc)
jacobian = fem.form(J)

A = fem.petsc.create_matrix(jacobian)
L = fem.petsc.create_vector(residual)

solver.setOperators(A)
dX = fem.Function(V)

# Newton iterations
i = 0
coords = V.tabulate_dof_coordinates()[:, 0]
sort_order = np.argsort(coords)
max_iterations = 25
solutions = np.zeros((max_iterations + 1, len(coords)))
solutions[0] = prc.x.array[sort_order]

# Solve problem in time steps
simu_time=simu_start_time+0.0001
while simu_time < end_time:  # Time steps
# Newton iterations
i = 0
converged=True
while i < max_iterations:
# Assemble Jacobian and residual
with L.localForm() as loc_L:
loc_L.set(0)
A.zeroEntries()
fem.petsc.assemble_matrix(A, jacobian,bcs=bc)
A.assemble()
fem.petsc.assemble_vector(L, residual)

# Scale residual by -1
L.scale(-1)

# Compute b - J(u_D-u_(i-1))
fem.petsc.apply_lifting(L, [jacobian], [bc], x0=[prc.vector], scale=1)
# Set dx|_bc = u_{i-1}-u_D
fem.petsc.set_bc(L, bc, prc.vector, 1.0)

# Solve linear problem
solver.solve(L, dX.vector)
dX.x.scatter_forward()
prc.x.array[:] += dX.x.array
i += 1

# Compute norm of update
correction_norm = dX.vector.norm(0)
# print(f"Iteration {i}: Correction norm {correction_norm}")
if correction_norm < 1e-10:
break
solutions[i,:] = prc.x.array[sort_order]
fem.petsc.assemble_vector(L, residual)
# print(f"Final residual {L.norm(0)}",flush=True)

# n, converged = solver.solve(prc)
prc.x.scatter_forward()
prc_prev.x.array[:]=prc.x.array
prc_prev.x.scatter_forward()
# if(simu_time%0.2<0.0001):
#     xdmf.write_function(prc, simu_time)
xdmf.write_function(prc, simu_time)
dt.value=set_time_step(dt.value,converged,i)
# print(dt.value,flush=True)
simu_time += dt.value

xdmf.close()
endT = time.time()
print("The time of execution of above program is :",(endT-startT),flush=True)
``````

Looked at this as well Point Sources (Redux) ?

1 Like

Thank you for this suggestion. I have looked into that post which should help.
However I think I am entangled into dimension inconsistencies as my attempt at integrating the provided code into my own throws ‘TypeError: determine_point_ownership(): incompatible function arguments’ when calling compute_cell_contributions.

MWE:

``````L = 1
W = 0.2
diffusion_time_step=0.01
end_time=1
simu_start_time=0
max_dif_tstep=0.2

import numpy as np
import ufl

from mpi4py import MPI
from petsc4py.PETSc import ScalarType

from dolfinx import mesh, fem, plot, io, nls, cpp
from dolfinx.fem import FunctionSpace,VectorFunctionSpace, Function,TensorFunctionSpace, Constant
from ufl import TestFunction, TrialFunction, dx, inner,grad, Measure,derivative,system
from petsc4py import PETSc
import time

startT = time.time()

domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0,0,0]), np.array([L, W, W])],
[30,30,30], cell_type=mesh.CellType.hexahedron)

V=VectorFunctionSpace(domain, ("CG", 1), dim=3)
V2=V.sub(0).collapse()[0]

p = TestFunction(V)
prc_prev = Function(V)
prc = Function(V)

dt=Constant(domain, diffusion_time_step)

def set_time_step(dt,converged,n):
if converged and n<5:
new_dt=min(max_dif_tstep,2*dt)
else: new_dt=dt
return new_dt

D_p = Constant(domain,2.4e3) # Diffusion coefficient for PTHrP
D_pr = Constant(domain,2.4e-15) # Diffusion coefficient for PTHrP
kap = Constant(domain,1e5) # Binding association rate for PTHrP-receptor
kdp = Constant(domain,1.86e-3)  # Disassociation rate for PTHrP-receptor
psi_p = Constant(domain,0.01) # Decay rate for PTHrP
psi_pr = Constant(domain,1e-50) # Decay rate for PTHrP receptor
psi_cp = Constant(domain,1e-15) # Decay rate for PTHrP-receptor Complex

xdmf = io.XDMFFile(domain.comm, "MiniExampleCustN.xdmf", "w")
xdmf.write_mesh(domain)

def compute_cell_contributions(V, points):
# Determine what process owns a point and what cells it lies within
mesh = V.mesh
_, _, owning_points, cells = 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 = fem.Expression(u, ref_x, comm=MPI.COMM_SELF)
values = expr.eval(mesh, np.asarray(cells, dtype=np.int32))

# 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=PETSc.ScalarType)
return cells, basis_values

if domain.comm.rank == 0:
points = np.array([[0.68, 0.36, 0]], dtype=domain.geometry.x.dtype) # ATTENTION here go the cell locations
else:
points = np.zeros((0, 3), dtype=domain.geometry.x.dtype)

cells, basis_values = compute_cell_contributions(V2, points)

# Supply term
#def sup_func(x):
#    return np.exp(-100*(x[0]%0.1+x[1]%0.1+x[2]%0.1))
p_i = Function(V2)
#p_i.interpolate(sup_func)
p_i.x.array[:] =0
for cell, basis_value in zip(cells, basis_values):
dofs = V2.dofmap.cell_dofs(cell)
p_i.x.array[dofs] += basis_value
p_sup = p_i

# PTHrP receptor presentation
def rec_func(x):
return np.exp(-10*(x[0]%0.15+x[1]%0.15+x[2]%0.15))
rp_i = Function(V.sub(0).collapse()[0])
rp_i.interpolate(rec_func)

subspace,submap=V.sub(1).collapse()
prc.x.array[submap]=rp_i.x.array
prc_prev.x.array[submap]=rp_i.x.array

# Variational problem
F = (
(- (prc[0] - prc_prev[0])*p[0]
- (prc[1] - prc_prev[1])*p[1]
- (prc[2] - prc_prev[2])*p[2])*dx()
+ (+ dt*p_sup*p[0])*dx()
+ (+ dt*(- kap*prc_prev[0]*prc[1]*p[0]
- kap*prc[0]*prc_prev[1]*p[1]
+ kap*prc_prev[0]*prc_prev[1]*p[2]
+ kdp*prc_prev[2]*p[0]
+ kdp*prc_prev[2]*p[1]
- kdp*prc[2]*p[2]))*dx()
+ (+ dt*(- psi_p*prc[0]*p[0]
- psi_pr*prc[1]*p[1]
- psi_cp*prc[2]*p[2]))*dx()
)

# Create boundary condition
tdim = V.mesh.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(V.mesh.topology)
fixed_dofs_0 = fem.locate_dofs_topological(V.sub(0), fdim, boundary_facets)
fixed_dofs_1 = fem.locate_dofs_topological(V.sub(1), fdim, boundary_facets)
fixed_dofs_2 = fem.locate_dofs_topological(V.sub(2), fdim, boundary_facets)
bc0 = fem.dirichletbc(PETSc.ScalarType(0), fixed_dofs_0, V.sub(0))
bc1 = fem.dirichletbc(PETSc.ScalarType(0), fixed_dofs_1, V.sub(1))
bc2 = fem.dirichletbc(PETSc.ScalarType(0), fixed_dofs_2, V.sub(2))
bc=[bc0,bc1,bc2]

# Define problem and custom Newton solver

solver = PETSc.KSP().create(domain.comm)
opts = PETSc.Options()
option_prefix =solver.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "hypre"
# opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
solver.setFromOptions()

residual = fem.form(F)
J = ufl.derivative(F, prc)
jacobian = fem.form(J)

A = fem.petsc.create_matrix(jacobian)
L = fem.petsc.create_vector(residual)

solver.setOperators(A)
dX = fem.Function(V)

# Newton iterations
i = 0
coords = V.tabulate_dof_coordinates()[:, 0]
sort_order = np.argsort(coords)
max_iterations = 25
solutions = np.zeros((max_iterations + 1, len(coords)))
solutions[0] = prc.x.array[sort_order]

# Solve problem in time steps
simu_time=simu_start_time+0.0001
while simu_time < end_time:  # Time steps
# Newton iterations
i = 0
converged=True
while i < max_iterations:
# Assemble Jacobian and residual
with L.localForm() as loc_L:
loc_L.set(0)
A.zeroEntries()
fem.petsc.assemble_matrix(A, jacobian,bcs=bc)
A.assemble()
fem.petsc.assemble_vector(L, residual)

# Scale residual by -1
L.scale(-1)

# Compute b - J(u_D-u_(i-1))
fem.petsc.apply_lifting(L, [jacobian], [bc], x0=[prc.vector], scale=1)
# Set dx|_bc = u_{i-1}-u_D
fem.petsc.set_bc(L, bc, prc.vector, 1.0)

# Solve linear problem
solver.solve(L, dX.vector)
dX.x.scatter_forward()
prc.x.array[:] += dX.x.array
i += 1

# Compute norm of update
correction_norm = dX.vector.norm(0)
# print(f"Iteration {i}: Correction norm {correction_norm}")
if correction_norm < 1e-10:
break
solutions[i,:] = prc.x.array[sort_order]
fem.petsc.assemble_vector(L, residual)
# print(f"Final residual {L.norm(0)}",flush=True)

# n, converged = solver.solve(prc)
prc.x.scatter_forward()
prc_prev.x.array[:]=prc.x.array
prc_prev.x.scatter_forward()
# if(simu_time%0.2<0.0001):
#     xdmf.write_function(prc, simu_time)
xdmf.write_function(prc, simu_time)
dt.value=set_time_step(dt.value,converged,i)
# print(dt.value,flush=True)
simu_time += dt.value

xdmf.close()
endT = time.time()
print("The time of execution of above program is :",(endT-startT),flush=True)

``````

For future posts, please post what version of DOLFINx that you are running. As far as I can tell, you are running v0.6.0, where you need to adapt the code to the old API:

``````L = 1
W = 0.2
diffusion_time_step=0.01
end_time=1
simu_start_time=0
max_dif_tstep=0.2

import numpy as np
import ufl

from mpi4py import MPI
from petsc4py.PETSc import ScalarType

from dolfinx import mesh, fem, plot, io, nls, cpp
from dolfinx.fem import FunctionSpace,VectorFunctionSpace, Function,TensorFunctionSpace, Constant
from ufl import TestFunction, TrialFunction, dx, inner,grad, Measure,derivative,system
from petsc4py import PETSc
import time

startT = time.time()

domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0,0,0]), np.array([L, W, W])],
[30,30,30], cell_type=mesh.CellType.hexahedron)

V=VectorFunctionSpace(domain, ("CG", 1), dim=3)
V2=V.sub(0).collapse()[0]

p = TestFunction(V)
prc_prev = Function(V)
prc = Function(V)

dt=Constant(domain, diffusion_time_step)

def set_time_step(dt,converged,n):
if converged and n<5:
new_dt=min(max_dif_tstep,2*dt)
else: new_dt=dt
return new_dt

D_p = Constant(domain,2.4e3) # Diffusion coefficient for PTHrP
D_pr = Constant(domain,2.4e-15) # Diffusion coefficient for PTHrP
kap = Constant(domain,1e5) # Binding association rate for PTHrP-receptor
kdp = Constant(domain,1.86e-3)  # Disassociation rate for PTHrP-receptor
psi_p = Constant(domain,0.01) # Decay rate for PTHrP
psi_pr = Constant(domain,1e-50) # Decay rate for PTHrP receptor
psi_cp = Constant(domain,1e-15) # Decay rate for PTHrP-receptor Complex

xdmf = io.XDMFFile(domain.comm, "MiniExampleCustN.xdmf", "w")
xdmf.write_mesh(domain)

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

# Pull owning points back to reference cell
mesh_nodes = domain.geometry.x
cmap = domain.geometry.cmap
ref_x = np.zeros((len(cells), domain.geometry.dim),
dtype=domain.geometry.x.dtype)
for i, (point, cell) in enumerate(zip(owning_points, cells)):
geom_dofs = domain.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 = fem.Expression(u, ref_x, comm=MPI.COMM_SELF)
values = expr.eval(domain, np.asarray(cells, dtype=np.int32))

# 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=PETSc.ScalarType)
return cells, basis_values

if domain.comm.rank == 0:
points = np.array([[0.68, 0.36, 0]], dtype=domain.geometry.x.dtype) # ATTENTION here go the cell locations
else:
points = np.zeros((0, 3), dtype=domain.geometry.x.dtype)

cells, basis_values = compute_cell_contributions(V2, points)

# Supply term
#def sup_func(x):
#    return np.exp(-100*(x[0]%0.1+x[1]%0.1+x[2]%0.1))
p_i = Function(V2)
#p_i.interpolate(sup_func)
p_i.x.array[:] =0
for cell, basis_value in zip(cells, basis_values):
dofs = V2.dofmap.cell_dofs(cell)
p_i.x.array[dofs] += basis_value
p_sup = p_i

# PTHrP receptor presentation
def rec_func(x):
return np.exp(-10*(x[0]%0.15+x[1]%0.15+x[2]%0.15))
rp_i = Function(V.sub(0).collapse()[0])
rp_i.interpolate(rec_func)

subspace,submap=V.sub(1).collapse()
prc.x.array[submap]=rp_i.x.array
prc_prev.x.array[submap]=rp_i.x.array

# Variational problem
F = (
(- (prc[0] - prc_prev[0])*p[0]
- (prc[1] - prc_prev[1])*p[1]
- (prc[2] - prc_prev[2])*p[2])*dx()
+ (+ dt*p_sup*p[0])*dx()
+ (+ dt*(- kap*prc_prev[0]*prc[1]*p[0]
- kap*prc[0]*prc_prev[1]*p[1]
+ kap*prc_prev[0]*prc_prev[1]*p[2]
+ kdp*prc_prev[2]*p[0]
+ kdp*prc_prev[2]*p[1]
- kdp*prc[2]*p[2]))*dx()
+ (+ dt*(- psi_p*prc[0]*p[0]
- psi_pr*prc[1]*p[1]
- psi_cp*prc[2]*p[2]))*dx()
)

# Create boundary condition
tdim = V.mesh.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(V.mesh.topology)
fixed_dofs_0 = fem.locate_dofs_topological(V.sub(0), fdim, boundary_facets)
fixed_dofs_1 = fem.locate_dofs_topological(V.sub(1), fdim, boundary_facets)
fixed_dofs_2 = fem.locate_dofs_topological(V.sub(2), fdim, boundary_facets)
bc0 = fem.dirichletbc(PETSc.ScalarType(0), fixed_dofs_0, V.sub(0))
bc1 = fem.dirichletbc(PETSc.ScalarType(0), fixed_dofs_1, V.sub(1))
bc2 = fem.dirichletbc(PETSc.ScalarType(0), fixed_dofs_2, V.sub(2))
bc=[bc0,bc1,bc2]

# Define problem and custom Newton solver

solver = PETSc.KSP().create(domain.comm)
opts = PETSc.Options()
option_prefix =solver.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "hypre"
# opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
solver.setFromOptions()

residual = fem.form(F)
J = ufl.derivative(F, prc)
jacobian = fem.form(J)

A = fem.petsc.create_matrix(jacobian)
L = fem.petsc.create_vector(residual)

solver.setOperators(A)
dX = fem.Function(V)

# Newton iterations
i = 0
coords = V.tabulate_dof_coordinates()[:, 0]
sort_order = np.argsort(coords)
max_iterations = 25
solutions = np.zeros((max_iterations + 1, len(coords)))
solutions[0] = prc.x.array[sort_order]

# Solve problem in time steps
simu_time=simu_start_time+0.0001
while simu_time < end_time:  # Time steps
# Newton iterations
i = 0
converged=True
while i < max_iterations:
# Assemble Jacobian and residual
with L.localForm() as loc_L:
loc_L.set(0)
A.zeroEntries()
fem.petsc.assemble_matrix(A, jacobian,bcs=bc)
A.assemble()
fem.petsc.assemble_vector(L, residual)

# Scale residual by -1
L.scale(-1)

# Compute b - J(u_D-u_(i-1))
fem.petsc.apply_lifting(L, [jacobian], [bc], x0=[prc.vector], scale=1)
# Set dx|_bc = u_{i-1}-u_D
fem.petsc.set_bc(L, bc, prc.vector, 1.0)

# Solve linear problem
solver.solve(L, dX.vector)
dX.x.scatter_forward()
prc.x.array[:] += dX.x.array
i += 1

# Compute norm of update
correction_norm = dX.vector.norm(0)
# print(f"Iteration {i}: Correction norm {correction_norm}")
if correction_norm < 1e-10:
break
solutions[i,:] = prc.x.array[sort_order]
fem.petsc.assemble_vector(L, residual)
# print(f"Final residual {L.norm(0)}",flush=True)

# n, converged = solver.solve(prc)
prc.x.scatter_forward()
prc_prev.x.array[:]=prc.x.array
prc_prev.x.scatter_forward()
# if(simu_time%0.2<0.0001):
#     xdmf.write_function(prc, simu_time)
xdmf.write_function(prc, simu_time)
dt.value=set_time_step(dt.value,converged,i)
# print(dt.value,flush=True)
simu_time += dt.value

xdmf.close()
endT = time.time()
print("The time of execution of above program is :",(endT-startT),flush=True)
``````

Please note that `determine_point_ownership` has been vastly improved in v0.7.x