Hi everyone,
I have a coupled heat equation here It works for some mesh, as in
# Mesh - time
T = 15 # final time
num_steps = 50 # number of time steps
dt = T/num_steps # time step size
# Mesh - space ---- 1D -> 2D
n_x = 200 # number of pts
n_y = 200 # number of pts
it runs and produces a result.
However when I change the above parameters as
# Mesh - time
T = 150 # final time
num_steps = 500 # number of time steps
dt = T/num_steps # time step size
# Mesh - space ---- 1D -> 2D
n_x = 20 # number of pts
n_y = 20 # number of pts
It gives the following error:
Traceback (most recent call last):
File "trwd2S.py", line 89, in <module>
solve(F == 0, u, BCs)
File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py", line 220, in solve
_solve_varproblem(*args, **kwargs)
File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py", line 266, in _solve_varproblem
solver.solve()
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to solve nonlinear system with NewtonSolver.
*** Reason: Newton solver did not converge because maximum number of iterations reached.
*** Where: This error was encountered inside NewtonSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset: 74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------
# nonlinear 2 pde system
# Packages
from dolfin import *
from mshr import *
import numpy as np
# Mesh - time
T = 150 # final time
num_steps = 500 # number of time steps
dt = T/num_steps # time step size
# Mesh - space ---- 1D -> 2D
n_x = 20 # number of pts
n_y = 20 # number of pts
x_min = -50.0 # min position in x
x_max = +50.0 # max position in x
y_min = -50.0 # min position in y
y_max = +50.0 # max position in y
mesh = RectangleMesh(Point(-50, -50), Point(+50, +50), n_x, n_y, "right")
# Finite elements - 1D -> 2D triangle
P1 = FiniteElement('CG',triangle, 1)
V = FunctionSpace(mesh, MixedElement([P1, P1]))
# Unknown Functions - set & split
v = TestFunction(V)
u = Function(V) # since it is a nonlinear function -> no need for TrialFunction just Function can be used. WHY?
v1, v2 = split(v)
u1, u2 = split(u)
# Initial condition ---- 1D -> 2D +pow(x[1],2))
u_0 = Expression(('0.20017 + 2*exp(-50*(pow(x[0],2)+pow(x[1],2)))','0.961477'), degree=2) # WHY DEGREE IS TWO?
# Known Function - set & split (from interpolating IC)
u_n = interpolate(u_0,V)
u_n1, u_n2 = split(u_n)
# Boundary - position ---- 1D -> 2D
def boundary(x):
return near(x[0], x_min, DOLFIN_EPS) or near(x[0], x_max, DOLFIN_EPS) or near(x[1], x_min, DOLFIN_EPS) or near(x[1], x_max, DOLFIN_EPS)
# Boundary condition
u10 = Expression('0.20017', degree=2)
u20 = Expression('0.961477', degree=2)
bc1 = DirichletBC(V.sub(0), u10, boundary)
bc2 = DirichletBC(V.sub(1), u20, boundary)
BCs = [bc1, bc2]
# Equation - coefficients
b = Constant(0.111)
mu = Constant(0.28900)
D = Constant(0.00001)
Dc = Constant(20.0)
D0 = Constant(0.1)
T = Constant(5.71429)
K = Constant(0.142857)
K1 = Constant(46.2857)
K2 =Constant(1.0)
# Equation - Crank-Nicolson
F1 = u1*v1*dx - u_n1*v1*dx + 0.5*(dt*D0*dot(grad(u1), grad(v1))*dx - dt*v1*(mu*K1*((b+(u1))/(1+(u1)))*u2 - ((T*(u1))/(K+(u1))))*dx
+ dt*D0*dot(grad(u_n1), grad(v1))*dx- dt*v1*(mu*K1*((b+(u_n1))/(1+(u_n1)))*u_n2 - ((T*(u_n1))/(K+(u_n1))))*dx)
F2 = u2*v2*dx - u_n2*v2*dx + 0.5*(dt*D*dot(grad(u2), grad(v2))*dx - dt*((1/(1+pow(u1,2)) - u2))*v2*dx
+ dt*D*dot(grad(u_n2), grad(v2))*dx - dt*v2*((1/(1+pow(u_n1,2)) - u_n2))*dx)
F = F1 + F2
# Create Solution Files
file_u1 = File("nonlin/var1_nonlin.pvd")
file_u2 = File("nonlin/var2_nonlin.pvd")
# Time - created
t = 0.0
# Time - loop
for n in range(num_steps):
# Update current time
t += dt
# Solve PDE system at time step
solve(F == 0, u, BCs)
# Save solution at time step
_u_1, _u_2 = u.split()
file_u1 << (_u_1, t)
file_u2 << (_u_2, t)
# Reassign previous time step
u_n.assign(u)
I do somethings without understanding the reason behind it, I write below incase they might be creating a problem.
FiniteElement('CG',triangle, 1) ( ‘CG’&1.It has something to do with differential forms thats all I know.)
u10 = Expression('0.20017', degree=2) (why the degree is 2 ? Honestly have no idea.)
Crank-Nicolson method (I only use this someone told me so)
Bonus question: Can you suggest a source for learning these? I wish to implement these without too much delving into theory. (I lost in pages of books, maybe someone have a compact lecture note, or an introduction of a thesis where these things are summarized.)
Thank you all in advance!