### Point source for nonlinear problem

315

views

1

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.

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

### 2 Answers

3

Thanks to MiroK :

https://fenicsproject.org/qa/2893/point-source-in-rhs-of-heat-eq

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

https://fenicsproject.org/qa/2893/point-source-in-rhs-of-heat-eq

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:
bc.apply(u.vector())
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)
bc.homogenize()
bcs_u.append(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 :
P.apply(b)
# 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.1

This seems to work. Thank you!

written
7 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
7 months ago by
Thomas Roy

Whoops! forgot to apply the boundary condition to "u" before entering Newton... See updated code above.

written
7 months ago by
pf4d

1

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
else:
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=Wells(degree=3)
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
7 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
7 months ago by
Ruben Gonzalez

No, what I said was that you

To do this you can use:

https://fenicsproject.org/olddocs/dolfin/1.6.0/python/programmers-reference/cpp/fem/FiniteElement.html#dolfin.cpp.fem.FiniteElement.evaluate_basis

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.

*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:

https://fenicsproject.org/olddocs/dolfin/1.6.0/python/programmers-reference/cpp/fem/FiniteElement.html#dolfin.cpp.fem.FiniteElement.evaluate_basis

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
7 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
7 months ago by
pf4d

Please login to add an answer/comment or follow this question.

Look at Farrell et al. paper for further details.