Hello, I’m new to FEniCS, please help. Why at different points in time do the graphs have tails that tend to “2”? Maybe it’s in the boundary conditions?
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
import scipy.special
import math
import os
plt.figure(figsize=(10, 5), dpi=80)
set_log_level(LogLevel.ERROR)
alpha = 1 #Показатель производной
T = 1 #Максимальное время (до которого рассчитываем u)
M = 100 #Количество интервалов во временной сетке
dt = T/M #Шаг по времени
nx = 65
t = 0
k0 = 1
A = 1
sigma = 2
n = -1/(0.5*sigma) ### n = -1/sigma - Sрежим, n < -1/sigma - HSрежим, 0 > n > -1/sigma - LSрежим
print(f'n={n}')
L_f = np.sqrt(2*k0*(A**sigma)*(sigma+2)/sigma)
mesh = IntervalMesh(nx, 0, L_f) #Задаем пространственную сетку
V = FunctionSpace(mesh, 'CG', 10) #Задаем пространство функций
xi = Expression('x[0]/(pow(k0, 0.5)*pow(A, sigma*0.5)*pow(1-t, (1+n*sigma)*0.5))', degree=10, sigma=sigma, n=n, t=t, k0=k0, A=A)
if n == -1/sigma:
u0 = Expression('A*pow(1-t,n)*pow((1-(x[0])/L_f),(2/sigma))', degree=10, sigma=sigma, t=t, n=n, A=A, k0=k0, L_f=L_f) #точное решение
else:
print('u02')
u0 = Expression('A*pow(1-t,n)*pow((1-(xi)/L_f),(2/sigma))', degree=10, sigma=sigma, t=t, n=n, A=A, k0=k0, xi=xi, L_f=L_f) #точное решение
#Задаем граничные условия
def b1(x, on_boundary):
return near(x[0], 0)
def b2(x, on_boundary):
if x[0] <= L_f:
return near(x[0], L_f)
else:
return 0
bc1 = DirichletBC(V, u0, b1)
bc2 = DirichletBC(V, 0, b2)
#Вводим функции u и v, для того чтобы в дальнейшем определить a(u,v) и L(v)
v = TestFunction(V)
u = Function(V)
#Функции q(u) и f(u) из постановки задачи
def q(uk):
return uk**sigma
ut = [interpolate(u0, V) for j in range(M)] #Значение функции u в начальный момент времени
for k in range(1, M):
t = dt * k
xi.t = t
u0.t = t
u0.xi = xi
F = u*v*dx + (dt**alpha)*dot(q(ut[k-1])*grad(u), grad(v))*dx - (ut[0] - sum( (-1)**j * scipy.special.binom(alpha, j) * (ut[k-j]-ut[0]) for j in range(1, k)) )*v*dx
J = derivative(F, u)
problem = NonlinearVariationalProblem(F, u, [bc1,bc2], J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 1E-5
prm['newton_solver']['relative_tolerance'] = 1E-4
prm['newton_solver']['maximum_iterations'] = 50000
prm['newton_solver']['relaxation_parameter'] = 1.0
solver.solve()
ut[k].assign(u)
if k % 10 == 0 or k == M - 1:
plot(u,label='t=%.2f'%(t))
plt.title(f'alpha={alpha}, \u0394x={1/nx}')
plt.grid(True)
plt.ylim(-0.5, 25)
plt.show()
print('complete')