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)