Hello all,

The code doesn’t give an error and it doesn’t evolve with time.

Below there is both a `PointSource`

and a velocity vector. Before I add PointSource, The code was running good. I am ran out of ideas.

There might be a problem due to solver, or some other thing that I couldn’t realize.

Also, where can I learn how to pick an appropriate solver for an equation? I look for “Numerical Linear Algebra with Applications” sources (books lecture notes etc.) but if you would suggest one that would be great?

What I did was the following:

```
#BUILD2 - pointsource
# PART - 1/2
# assemble before applying pointsource
A, b = assemble_system(a, L, bcs)
delta = PointSource(V, Point(4., 4.), -1)
delta.apply(b)
# PART - 2/2
# if you define solution as
uh = Function(V)
# Solve PDE system at time step
solve(A, uh.vector(), b)
```

before I was using

```
# Solve variational problem for time step
solve(a == L, u, bcs)
```

The whole code:

```
# Import Packages
from dolfin import *
from mshr import *
import time
start_time = time.time()
# Time Stepping
tmax = 60
num_steps = 200 # number of time steps
dt = tmax/num_steps # time step size
# meshs
mesh = RectangleMesh(Point(0.0, 0.0), Point(8.0, 8.0), 160, 160, "right") #1D -> 2D
P1 = FiniteElement('CG', triangle, 1)
# Function space & trial function - concentration
V = FunctionSpace(mesh, P1)
u = TrialFunction(V)
v = TestFunction(V)
# BUILD1- beta p1
# Define function space for velocity
W = VectorFunctionSpace(mesh, 'P', 2)
w = Function(W)
# velocity
beta = project(Expression(('-2','0'), degree=1), W)
# Initial condition #1D -> 2D
u_0 = Expression(('3.0'), degree=1)
u_n = interpolate(u_0, V)
# Boundary #1D -> 2D
def boundary(x):
return near(x[0], 0.0, DOLFIN_EPS) or near(x[0], 8.0, DOLFIN_EPS) or near(x[1], 0.0, DOLFIN_EPS) or near(x[1], 8.0, DOLFIN_EPS)
# Condition
u10 = Expression('1.0', degree=1)
bcs = DirichletBC(V, u10, boundary)
# PDE - Crank Nicolson
F = (u - u_n)/dt*v*dx + 0.5*(Constant(0.01)*dot(grad(u), grad(v))*dx + dot(beta, grad(u))*v*dx + Constant(2.0)*v*dx
+ Constant(0.01)*dot(grad(u_n), grad(v))*dx + dot(beta, grad(u_n))*v*dx + Constant(2.0)*v*dx)
a = lhs(F)
L = rhs(F)
# assemble before applying pointsource
A, b = assemble_system(a, L, bcs)
delta = PointSource(V, Point(4., 4.), +4)
delta.apply(b)
# Solution
u = Function(V)
# Create VTK files for visualization output
vtkfile = File('bbb/bbb.pvd')
# Time-stepping
t = 0
for n in range(num_steps):
# Update current time
t += dt
# Solve variational problem for time step
solve(A, u.vector(), b, 'bicgstab', 'petsc_amg')
# Save solution to file (VTK)
vtkfile << (u, t)
# Update previous solution
u_n.assign(u)
# Print update
print ("\nTime " + str(t) + " of " + str(tmax) + "\n")
```

Any suggestion would be appreciated,

Bests.