Code doesn't evolve with time - PointSource related

1

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.