How to correctly apply a point source in nonlinear problem in fenicsx

dolfinx
1

I am simulating a problem where flow injected through a point source on a horizontal plane. The axisymmetric dimensionless equation is

\begin{aligned} xu_t+\left(x\partial_x+1\right)\Gamma&=\frac{\delta(x)}{2\pi } \quad\text{in } \Omega=[0, 1]\\ u(0,x)&=0\\ u_x(t,0)&=0,\quad u(t,1)=0 \end{aligned}

where the flux \Gamma=-u^3u_x and the weak form is

\int_0^1 \left(vxu_t-x\Gamma v_x\right)\mathrm{d} x=\frac{v(0)}{2\pi}

I checked this post on point source, but that is a linear problem and the v(0)/2\pi term is added to the RHS b. Is there a way to add a point source for a nonlinear problem? The following MWE is for a simplified problem

\begin{aligned} xu_t+\left(x\partial_x+1\right)\Gamma&=0\quad\text{in } \Omega=[0, 1]\\ u(0,x)&=\frac{1}{1+100x^2}\\ u_x(t,0)&=0,\quad u(t,1)=0 \end{aligned} 
"""
Time-dependent 1D lubrication equation [@garel2014analogue]
u_t+Γ_x+Γ/x=0, Γ_x=-u^3 u_x
u_x(0,t)=0, u(R,t)=0, u(0,t)=exp(-100 x^2)

The weak formulation is
∫(vxu_t-xΓv_x)dx = 0
"""

# %%
import matplotlib.pylab as plt
import numpy as np
from dolfinx import mesh
from dolfinx.fem import (
    Constant,
    Function,
    dirichletbc,
    form,
    functionspace,
    locate_dofs_geometrical,
)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.io import XDMFFile
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from ufl import TestFunction, derivative, dx, grad, inner


# Define temporal parameters
T = 1  # final time
num_steps = 100  # time steps
dt = T / num_steps

# Create mesh and define function space
nx = 200  # number of intervals in space
rank = MPI.COMM_WORLD.Get_rank()
domain = mesh.create_interval(MPI.COMM_WORLD, nx, [0, 1])
V = functionspace(domain, ("Lagrange", 1))

# Create initial condition
initial_condition = lambda x: 1 / (1 + 100 * x[0] ** 2)
# Define solution variable, and interpolate initial solution for visualization in Paraview
u = Function(V)
u_n = Function(V)
u.interpolate(initial_condition)
u.x.scatter_forward()
u_n.x.array[:] = u.x.array


# Define boundary conditions
bc_right = dirichletbc(
    Constant(domain, PETSc.ScalarType(0)),
    locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 1)),
    V,
)
bcs = [bc_right]

# Time-dependent output
# xdmf = XDMFFile(domain.comm, "result/eqn02_diffusion.xdmf", "w")
# xdmf.write_mesh(domain)

# xdmf.write_function(uh, t)

v = TestFunction(V)

F = inner(u - u_n, v) * dx + dt * u**3 * inner(grad(u), grad(v)) * dx


# Compute the Jacobian (derivative of F with respect to u)
J = derivative(F, u)

problem = NonlinearProblem(F, u, bcs=bcs, J=J)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-8

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()

x = V.tabulate_dof_coordinates()[:, 0]

# Updating the solution and right hand side per time step
t = 0
for i in range(num_steps):
    t += dt

    # Solve linear problem
    n, converged = solver.solve(u)
    if rank == 0:
        print(f"Number of interations: {n:d}")
    u.x.scatter_forward()
    u_n.x.array[:] = u.x.array
    u.x.scatter_forward()

    # Update solution at previous time step (u_n)
    u_n.x.array[:] = u.x.array

    # Write solution to file
    # xdmf.write_function(uh, t)
    if i % 10 == 0:
        plt.plot(x, u.x.array)
plt.show()
# xdmf.close()
# %%