Hey,
wanted to use the function circle = Circle (0,0,1) and mesh = Mesh(circle,20) to solve the simple parabolic PDE $\partial_t u - \Delta u = 0$ with initial condition $u(0,x_1, x_2) = x_1 + 3$, but the NewtonSolver is not converging.
However, for mesh = UnitSquareMesh(20, 20) or mesh = UnitCircle(20) everything is fine and I get the expected steady state solution $u = 3.5$, respectively $u = 3.0$. What is going on, do I make some stupid mistake while using Circle()?
Thanks in advance for your help.
from __future__ import division
from dolfin import *
import numpy as np
import matplotlib as mpl
mpl.use( "agg" )
import matplotlib.pyplot as plt
from matplotlib import animation
from math import exp
from math import sin
from math import cos
from math import pow
from math import tanh
# to shut off the output
set_log_level(WARNING)
nx = 20
circle = Circle (0.0,0.0,1)
mesh = Mesh(circle,20)
#mesh = UnitSquareMesh(nx, nx)
#mesh = UnitCircle(20)
plot(mesh,axes=True)
interactive()
# definition of solving method
V = FunctionSpace(mesh, 'CG', 1)
T = 10
dt = 1/(80)
# ------------------- RHS, exact solution
# RHS
def rhs(u):
return 0
u_start = Expression("x[0] + 3")
# ------------------- RHS, exact solution
# ------------------- boundary conditions
# Define boundary segments
class outer_boundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-13 # tolerance for coordinate comparisons
return on_boundary and abs( x[0]*x[0] + x[1]*x[1] - 9) < tol
class left(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-13 # tolerance for coordinate comparisons
return on_boundary and abs(x[0]) < tol
z_old = Function(V)
z_old = interpolate(u_start,V)
z = Function(V)
z = interpolate(u_start,V)
dz = TrialFunction(V)
w = TestFunction(V)
# - dt*g_N*w*dssE(2)
F_N = inner(z-z_old,w)*dx + dt*dot(grad(z), grad(w))*dx - inner(dt*rhs(z),w)*dx
dF_N = derivative(F_N, z, dz)
problem = NonlinearVariationalProblem(F_N, z, J = dF_N)
pdesys_Euler = NonlinearVariationalSolver(problem)
plot(z, title='impl Euler',axes=True)
interactive()
# ---------- update of time dependent data
t = dt
# ---------- time steps
while t <= T+dt:
# ---------- Newton
# Euler
pdesys_Euler.solve()
# ---------- assigning of solution from previous time step
z_old.assign(z)
if (t >0.99 and t < 1.01):
plot(z)
interactive()
if (t >1.99 and t < 2.01):
plot(z)
interactive()
t += dt
