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[0], x_min, tol)
class BoundaryYMIN(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x[1], x_min, tol)
class BoundaryXMAX(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x[0], x_max, tol)
class BoundaryYMAX(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x[1], 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[0], betaexp[1])) # 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