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*inner(D_p*grad(prc[0]), grad(p[0])))*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)
                L.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
                
                # Scale residual by -1
                L.scale(-1)
                L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)

                # 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)
                L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)

                # 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*inner(D_p*grad(prc[0]), grad(p[0])))*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)
                L.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
                
                # Scale residual by -1
                L.scale(-1)
                L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)

                # 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)
                L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)

                # 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*inner(D_p*grad(prc[0]), grad(p[0])))*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)
                L.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
                
                # Scale residual by -1
                L.scale(-1)
                L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)

                # 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)
                L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)

                # 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

1 Like

Thank you for your reply.

I am actually running Dolfinx 0.7.0.0 (according to dolfinx.version command).
The code you provided still throws the same error for me:

TypeError: determine_point_ownership(): incompatible function arguments. The following argument types are supported:
1. (arg0: dolfinx.cpp.mesh.Mesh, arg1: numpy.ndarray[numpy.float64]) → Tuple[List[int], List[int], List[float], List[int]]

Thanks

Hello,
Sorry to come back to this question but I’m still struggling on that front… any idea how to address the issue?
Many thanks!

Upgrade to 0.7.2 or 0.7.3. They have critical bugfixes in that determine_point_ownership, see Release v0.7.2 · FEniCS/dolfinx · GitHub

Thank you. I upgraded to v0.8.0 and this now seems to work well.