 # PointSource with Robin BC - not calling Newton Solver

## Explanation

This is the first time I use Robin bc and PointSource.

The below code runs without an error yet, it doesn’t call Newton Solver. (I mean the Newton iterations that’s been counting on the terminal window. ) No error means no feedback, so I am stuck.

Any help would be deeply appreciated.

## Model

Below, a 2D model for advection diffusion reaction model

\partial_{t}{u} = 0.1\; \nabla^2 u - β . (\nabla u) + 0.9 \;δ(x-1.5) δ(y-1.5) - 0.5\; u

here \beta is a vector with β_{x} = 1.5 and β_{y} = 0.

The problem has a square domain (0,10)x(0,10) for t=(0,10).

Initial condition is
u (0,x,y) = 0.0

Boundary conditions are
\partial_{x} u(t,0,y)= β_{x} u(t,0,y)
\partial_{x} u(t,10,y)= β_{x} u(t,10,y)
\partial_{y} u(t,x,0)= β_{y} u(t,x,0)=0.0
\partial_{y} u(t,x,10)= β_{y} u(t,x,10)=0.0

Note: For the below BC definitions are taken as
Robin bc : -K \partial _{n} u = r (u - s)
Neum bc : -K \partial _{n} u = g

## Code

# Import Packages
from dolfin import *
from mshr import *
import numpy as np
import time

start_time = time.time()

# Parameters for domain
tmax = 10
xmax = 10

# Time Stepping
num_steps = 100       # number of time steps
dt = tmax/num_steps     # time step size

# Meshing & finite element                                                                     #1D -> 2D
n_x = 200       # Number of x Steps on Grid
n_y = 200       # Number of y Steps on Grid
# grid extremes are dimensionless
x_min = 0       # Minimum x on Grid
x_max = +xmax   # Maximum x on Grid
mesh = RectangleMesh(Point(x_min, x_min), Point(x_max, x_max), n_x, n_y, "right")               #1D -> 2D
P1 = FiniteElement('CG', triangle, 1)

# Function space & trial function - concentration
V = FunctionSpace(mesh, P1)
u = TrialFunction(V)
# Test function
v = TestFunction(V)

# Initial condition - value                                                               #1D -> 2D
u_0 = Expression(('0.0'), degree=1)
# Interpolate expression on to function space as first time step
u_n = interpolate(u_0,V)

# Subdomains
class BoundaryXMIN(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x, x_min, tol)
class BoundaryYMIN(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x, x_min, tol)
class BoundaryXMAX(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x, x_max, tol)
class BoundaryYMAX(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x, x_max, tol)

# Boundary splitted
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)

# Boundary subdomains
sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_markers.set_all(9999)
# Left
bxmin = BoundaryXMIN()
bxmin.mark(boundary_markers, 0)
# Top
bymax = BoundaryYMAX()
bymax.mark(boundary_markers, 1)
# Right
bxmax = BoundaryXMAX()
bxmax.mark(boundary_markers, 2)
# Bottom
bymin = BoundaryYMIN()
bymin.mark(boundary_markers, 3)

# Boundary values
g = Constant(0.0) # Neumann Value
r = Constant(-1.5) # Robin Value r(u-s)
s = Constant(0.0) # Robin Value r(u-s)

# Boundary dictionary
boundary_conditions = {0: {'Robin': (r, s)},
1: {'Neumann': g},
2: {'Robin': (r, s)},
3: {'Neumann': g}
}

# Collect boundaries
bcs = []
for j in boundary_conditions:
if 'Dirichlet' in boundary_conditions[j]:
bc = DirichletBC(V, boundary_conditions[j]['Dirichlet'],
boundary_markers, j)
bcs.append(bc)
integrals_N = []
for i in boundary_conditions:
if 'Neumann' in boundary_conditions[i]:
g = boundary_conditions[i]['Neumann']
integrals_N.append(g*v*ds(i))
integrals_R = []
for i in boundary_conditions:
if 'Robin' in boundary_conditions[i]:
r, s = boundary_conditions[i]['Robin']
integrals_R.append(r*(u - s)*v*ds(i))

# Measure
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

# Parameters
betaexp = Expression(("1.5", "0.0"), domain=mesh, degree = 1)
beta = as_vector((betaexp, betaexp))  # caution! double paranthesis!

# Files
file_u = File('adrfin/u.pvd')

# PDE - Crank Nicolson
F = (u - u_n)/dt*v*dx + sum(integrals_N) + sum(integrals_R) + 0.5*(Constant(0.1)*dot(grad(u), grad(v))*dx + dot(beta, grad(u))*v*dx + Constant(0.5)*u*v*dx
+ Constant(0.1)*dot(grad(u_n), grad(v))*dx + dot(beta, grad(u_n))*v*dx + Constant(0.5)*u_n*v*dx)

a = lhs(F)
L = rhs(F)
A, b = assemble_system(a, L, bcs)
delta = PointSource(V, Point(1.5, 1.5), -0.9)
delta.apply(b)

uh = Function(V)

t = 0.0
for n in range(num_steps):
# Solve PDE system at time step
solve(A, uh.vector(), b)
file_u << (uh, t)
# Reassign previous time step
u_n.assign(uh)
# Print update
print ("\nTime " + str(t) + " of " + str(tmax) + "\n")
# Update current time
t = t + dt


Hi @Kaa , can you clarify why you expect the Newton Solver to be invoked?

The Newton Solver is used for nonlinear systems, but you’ve invoked the high level solve() wrapper command, which calls the LinearSolver (a direct solver for the linear algebra, not iterative in the sense of NewtonSolver)

1 Like

Hi @dparsons,

I didn’t know the Newton Solver is only for nonlinear problems, that clears a lot. Thank you. Since I don’t see any evolution from the initial conditions and no “newton solver …” I thought they might be relevant.

Hi again, @dparsons ok the problem is not the solver but nothing evolves with time. I changed, i guess, every variable but it is still the same. -couldn’t change “not calling the newton solver” title to “nothing evolves with time”-

You had the right intuition. You want to update the initial guess each iteration from the previous solution, and you do that with the NewtonSolver. The problem is that you don’t have a NewtonSolver (a NonlinearVariationalSolver for a NonlinearVariationalProblem). You can’t provide an initial guess to the default LinearSolvers.

I have tried many things but none of them are working.
I really don’t know what is the problem here? I applied bcs to assembled pieces, used different preconditioners for linear solvers. Tried to use a nonlinear solver too, but it didn’t accepted my

F = b - A


since b is a vector and A is a matrix. So I wrote them as F = b*dx - A*dx I thought this would make them a form. It simply didn’t.

Any ideas???

## Trial1

# Tried to create a linear solver
# MAX_ITERS = 10000
# solver = PETScKrylovSolver('gmres','none')
# solver.ksp().setGMRESRestart(MAX_ITERS)
# solver.parameters["monitor_convergence"] = True
# solver.parameters['relative_tolerance'] = 1e-6
# solver.parameters['maximum_iterations'] = MAX_ITERS


## Trial2

uh = Function(V)
# Tried to create a nonlinear solver
# duh = TrialFunction(V)
# print(type(lhs(F)))
# print(type(rhs(F)))
# print(type(A))
# print(type(b))
# F = b - A
# Jac = derivative(F, uh, duh)
# problem = NonlinearVariationalProblem(F, uh, bcs, Jac)
# solver  = NonLinearVariationalSolver(problem)


## Trial3

 # solve(A, uh.vector(), b, "cg", "hypre_amg")


## Trial4

solve(A, uh.vector(), b, "cg", "ilu")


Below I left trials as they are with # So that one can see how the implemented in the code.

## The code

# Import Packages
from dolfin import *
from mshr import *
import numpy as np
import time

start_time = time.time()

# Parameters for domain
tmax = 600
xmax = 10

# Time Stepping
num_steps = 4000       # number of time steps
dt = tmax/num_steps     # time step size

# Meshing & finite element                                                                     #1D -> 2D
n_x = 50       # Number of x Steps on Grid
n_y = 50       # Number of y Steps on Grid
# grid extremes are dimensionless
x_min = 0       # Minimum x on Grid
x_max = +xmax   # Maximum x on Grid
mesh = RectangleMesh(Point(x_min, x_min), Point(x_max, x_max), n_x, n_y, "right")               #1D -> 2D
P1 = FiniteElement('CG', triangle, 1)

# Function space & trial function - concentration
V = FunctionSpace(mesh, P1)
u = TrialFunction(V)
# Test function
v = TestFunction(V)

# Initial condition - value                                                               #1D -> 2D
u_0 = Expression(('0.0'), degree=1)
# Interpolate expression on to function space as first time step
u_n = interpolate(u_0,V)

# Subdomains
class BoundaryXMIN(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x, x_min, tol)
class BoundaryYMIN(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x, x_min, tol)
class BoundaryXMAX(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x, x_max, tol)
class BoundaryYMAX(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x, x_max, tol)

# Boundary splitted
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)

# Boundary subdomains
sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_markers.set_all(9999)
# Left
bxmin = BoundaryXMIN()
bxmin.mark(boundary_markers, 0)
# Top
bymax = BoundaryYMAX()
bymax.mark(boundary_markers, 1)
# Right
bxmax = BoundaryXMAX()
bxmax.mark(boundary_markers, 2)
# Bottom
bymin = BoundaryYMIN()
bymin.mark(boundary_markers, 3)

# Boundary values
g = Constant(0.0) # Neumann Value
r = Constant(+10.5) # Robin Value r(u-s)#########################################
s = Constant(0.0) # Robin Value r(u-s)

# Boundary dictionary
boundary_conditions = {0: {'Robin': (r, s)},
1: {'Neumann': g},
2: {'Robin': (r, s)},
3: {'Neumann': g}
}

# Collect boundaries
bcs = []
for j in boundary_conditions:
if 'Dirichlet' in boundary_conditions[j]:
bc = DirichletBC(V, boundary_conditions[j]['Dirichlet'],
boundary_markers, j)
bcs.append(bc)
integrals_N = []
for i in boundary_conditions:
if 'Neumann' in boundary_conditions[i]:
g = boundary_conditions[i]['Neumann']
integrals_N.append(g*v*ds(i))
integrals_R = []
for i in boundary_conditions:
if 'Robin' in boundary_conditions[i]:
r, s = boundary_conditions[i]['Robin']
integrals_R.append(r*(u - s)*v*ds(i))

# Measure
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

# Parameters
betaexp = Expression(("1.5", "0.0"), domain=mesh, degree = 1) ##################################################
beta = as_vector((betaexp, betaexp))  # caution! double paranthesis!

# Files
file_u = File('adrfin/u.pvd')

# PDE - Crank Nicolson
F = (u - u_n)/dt*v*dx + sum(integrals_N) + sum(integrals_R) + 0.5*(Constant(1.0)*dot(grad(u), grad(v))*dx + dot(beta, grad(u))*v*dx + Constant(20.5)*u*v*dx #0.5)*u*v*dx
+ Constant(1.0)*dot(grad(u_n), grad(v))*dx + dot(beta, grad(u_n))*v*dx + Constant(20.5)*u_n*v*dx) #0.5)*u_n*v*dx)

a = lhs(F)
L = rhs(F)
A, b = assemble_system(a, L, bcs)
delta = PointSource(V, Point(1.5, 1.5), 30.0) #-0.9)
delta.apply(b)

# Tried to create a linear solver
# MAX_ITERS = 10000
# solver = PETScKrylovSolver('gmres','none')
# solver.ksp().setGMRESRestart(MAX_ITERS)
# solver.parameters["monitor_convergence"] = True
# solver.parameters['relative_tolerance'] = 1e-6
# solver.parameters['maximum_iterations'] = MAX_ITERS

uh = Function(V)
# Tried to create a nonlinear solver
# duh = TrialFunction(V)
# print(type(lhs(F)))
# print(type(rhs(F)))
# print(type(A))
# print(type(b))
# F = b*dx - A*dx
# Jac = derivative(F, uh, duh)
# problem = NonlinearVariationalProblem(F, uh, bcs, Jac)
# solver  = NonLinearVariationalSolver(problem)

t = 0.0
for n in range(num_steps):
# Solve PDE system at time step
# Tried to create a linear solver inside time loop
# solve(A, uh.vector(), b, "cg", "hypre_amg")
solve(A, uh.vector(), b, "cg", "ilu")
file_u << (uh, t)
# Reassign previous time step
u_n.assign(uh)
# Print update
print ("\nTime " + str(t) + " of " + str(tmax) + "\n")
# Update current time
t = t + dt



There’s an example demos of the Nonlinear solver at
https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/undocumented/contact-vi-snes/demo_contact-vi-snes.py