Point source for nonlinear problem

9 months ago by
I was wondering if there was a simple way to implement a point source term for nonlinear problems. For linear problems, we can use PointSource() and apply it to the rhs after assembling the system. This cannot be done when you have a NonlinearVariationalProblem.

I then defined a delta function using the limit definition. Obviously, this function is not well approximated by CG P1 elements. Would this work for DG?

My next idea is to modify the linearized system to add the source term, but I don't know how to do that.

Has anyone ever implemented a point source for a nonlinear problem? Thank you.
Community: FEniCS Project
Have you tried using vertex_integral instead of point source?
Look at Farrell et al. paper for further details.
written 9 months ago by Eleonora Piersanti  

2 Answers

9 months ago by
Thanks to MiroK :


for use of PointSource(), I have the manual Newton's method solved with a point source:
from fenics import *

mesh = UnitSquareMesh(50,50)
Q    = FunctionSpace(mesh, 'CG', 1)

u    = Function(Q)
du   = TrialFunction(Q)
v    = TestFunction(Q)
p    = Point(0.75,0.75)
P    = PointSource(Q, p, 100.0)
f    = Constant(100.0)

def boundary(x, on_boundary):  return on_boundary

bc   = DirichletBC(Q, 50.0, boundary)
bcs  = [bc]

R    = + inner(grad(u), grad(v)) * dx \
       - f * v * dx

J    = derivative(R, u, du)

# apply the boundary condition to the vector once, 
# the search direction is set to zero on the boundary
# with homogenize() below, and therefore 'u' will not 
# change along the Dirichlet boundaries :
for bc in bcs:
atol             = 1e-7
rtol             = 1e-10
relaxation_param = 1.0
max_iter         = 25
method           = 'mumps'
preconditioner   = 'default'
converged        = False
lmbda            = relaxation_param   # relaxation parameter
nIter            = 0                  # number of iterations
# need to homogenize the boundary, as the residual is always zero over
# essential boundaries :
bcs_u = []
for bc in bcs:
  bc = DirichletBC(bc)

# the direction of decent :
d = Function(u.function_space()) 

while not converged and nIter < max_iter:

  # assemble system :
  A, b    = assemble_system(J, -R, bcs_u)

  # apply the point source :

  # determine step direction :
  solve(A, d.vector(), b, method, preconditioner)

  # calculate residual :
  residual  = b.norm('l2')

  # set initial residual : 
  if nIter == 0:
    residual_0 = residual

  # the relative residual :
  rel_res = residual / (residual_0 + DOLFIN_EPS)

  # check for convergence :
  converged = residual < atol or rel_res < rtol
  # move U down the gradient :
  u.vector()[:] += lmbda*d.vector()
  # increment counter :
  nIter += 1

  # print info to screen :
  if MPI.rank(mpi_comm_world()) == 0:
    string = "Newton iteration %d: r (abs) = %.3e (tol = %.3e) " \
             +"r (rel) = %.3e (tol = %.3e)"
    print string % (nIter, residual, atol, rel_res, rtol)

File('u.pvd') << u

edit: correctly applied the boundary condition to the initial guess to fix issue described by below comment.
This seems to work. Thank you!
written 9 months ago by Thomas Roy  
The homogenization of the BCs just seems to remove my Dirichlet BCs altogether. Is there a way to fix that?
written 9 months ago by Thomas Roy  
Whoops! forgot to apply the boundary condition to "u" before entering Newton...  See updated code above.
written 9 months ago by pf4d  
9 months ago by
In a problem there are two wells, one injecting and one producing, PointSource() was not working for me, maybe an error in the implementation, any way, The solution was to setup a class that extends Expression

tol = 3E-5
class Wells(Expression):
    def eval(self,value,x):
        Q = 1E-2
        x0 = -25
        y0 = 0
        x1 = 25
        y1 = 0
        r0 = sqrt((x[0] - x0)**2+(x[1] - y0)**2)
        r1 = sqrt((x[0] - x1) ** 2 + (x[1] - y1) ** 2)

        if r0 < tol:
            value[0] = Q
        elif r1 < tol:
            value[0] = -Q
            value[0] = 0​

The tolerance gets a little tricky because, most likely, there is not a node in the position you want, a little bit of debugging is needed to check which one is the closest node and then adjust the tolerance.

To use it in your formulation you just do this, where q is your TestFunction


F = (w**3/(12*mu))*inner(grad(q),grad(p))*dx-f*q*dx

Hope this helps!!

Using this method, you could linearly interpolate a point located within the interior of a cell to its nodes; even better would be to create a mesh with nodes located exactly at the point source location.  GMSH will do this for you, perhaps MSHR?
written 9 months ago by pf4d  
It does not, because you are choosing the closest node to your desired point and assign the value, all others will be 0. you can verify that with
f_temp = assemble(f*q*dx)​

f_temp will be a vector where only your chosen points are nonzero.

To add accuracy, a mesh was created with GMSH, then manually modify the value of the closest point to the desired value and then convert to Fenics.

written 9 months ago by Ruben Gonzalez  
No, what I said was that you could interpolate the value of a point source located within the cell to each of that cell's nodes, not that your method does this.

To do this you can use:


And regarding GMSH, you can just specify the node locations exactly, then refine the interior points around these nodes you have chosen.  Once the mesh has been generated by GMSH, you will have a .msh file which you convert directly to FEniCS with "dofin-convert".  Then, when you use your method above to set the closest node, you will specify the position of the point source to machine precision.

Finally, I think that MSHR may already provide this functionality, but I am too lazy to find out for you.
written 9 months ago by pf4d  
Also, If you only have one point source, you could just adjust the node spacing and domain coordinates such that a node of a structured mesh will exist exactly where you want it.
written 9 months ago by pf4d  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »