FEniCS: UMFPACK сообщает, что решаемая матрица является единственной
Я создал разделенную область с уравнением Стокса в первом поддомене и смешанным уравнением Пуассона (Дарси) во втором поддомене. Я работаю с UnitSquare и поддомен 1 должен иметь интервал от 0 до 0,5 и поддомен 2 от 0,5 до 1.
Но теперь я получаю следующую ошибку:
Решение линейной вариационной задачи. Проблема UMFPACK, связанная с вызовом числового значения * Предупреждение: UMFPACK сообщает, что решаемая матрица является единственной. Проблема UMFPACK, связанная с вызовом для решения * Предупреждение: UMFPACK сообщает, что решаемая матрица является единственной. assert vmax>=vmin, "пустой диапазон, укажите vmin и / или vmax" Ошибка подтверждения: пустой диапазон, укажите vmin и / или vmax
Кто-нибудь может помочь? Спасибо!
Вот код:
enter code here
#-*- coding: utf-8 -*-
from dolfin import *
import numpy as np
# Define mesh
mesh = UnitSquare(32,32)
#Subdomain 1
# Gitter übergeben
subdomains = CellFunction("uint", mesh)
# Klasse des Teilgebiets
class Domain_1(SubDomain):
def inside(self, x, on_boundary):
return between(x[0], (0, 0.5)) # Koordinatenangabe des Teilgebiets
# Objekt der Klasse erstellen
sub_domain1 = Domain_1()
sub_domain1.mark(subdomains,0)
# Definition Funktionenräume
U = FunctionSpace(mesh, "CG", 2)
V = FunctionSpace(mesh, "CG", 1)
W = U*V
# Definition Trial- und Testfunktion
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
# boundary condition
p_in = 1
p_out = 0
noslip = DirichletBC(W.sub(0), (0),
"on_boundary && \
(x[1] <= DOLFIN_EPS | x[1] >= 0.5-DOLFIN_EPS)")
inflow = DirichletBC(W.sub(1), p_in, "x[0] <= 0.0 + DOLFIN_EPS*1000")
outflow = DirichletBC(W.sub(1), p_out, "x[0] >= 0.5 - DOLFIN_EPS*1000")
bcp = [noslip,inflow, outflow]
# Definition f
f = Expression("0")
# Variationsformulierung
a = inner(grad(u), grad(v))*dx + div(v)*p*dx(0) + q*div(u)*dx(0)
L = inner(f,v)*dx(0)
# Lösung berechnen
w = Function(W)
problem = LinearVariationalProblem(a, L, w, bcp)
solver = LinearVariationalSolver(problem)
solver.solve()
(u, p) = w.split()
# Subdomain 2
# Gitter übergeben
subdomains = CellFunction("uint", mesh)
# Klasse des Teilgebiets
class Domain_2(SubDomain):
def inside(self,x,on_boundary):
return between(x[0], (0.5,1.0)) # Koordinatenangabe des Teilgebiets
# Objekt der Klasse erstellen
sub_domain2 = Domain_2()
sub_domain2.mark(subdomains,1)
# Define function spaces and mixed (product) space
BDM = FunctionSpace(mesh, "BDM", 1)
DG = FunctionSpace(mesh, "DG", 0)
CG = FunctionSpace(mesh, "CG", 1)
W = MixedFunctionSpace([BDM, DG, CG])
# Define trial and test functions
(sigma, u, p) = TrialFunctions(W)
(tau, v, q) = TestFunctions(W)
#Define pressure boundary condition
p_in = 1
p_out = 0
noslip = DirichletBC(W.sub(1), (0),
"on_boundary && \
(x[1] <= 0.5 + DOLFIN_EPS | x[1] >= 1.0-DOLFIN_EPS)")
inflow = DirichletBC(W.sub(2), p_in, "x[0] <= 0.5 + DOLFIN_EPS*1000")
outflow = DirichletBC(W.sub(2), p_out, "x[0] >= 1.0 - DOLFIN_EPS*1000")
bcp = [noslip,inflow, outflow]
# Define f
#f = Expression("0")
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
# Define variational form
a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx(1) + inner(p,q)*dx(1) + u*q*dx(1)
L = f*v*dx(1)
# Compute solution
w = Function(W)
problem = LinearVariationalProblem(a, L, w, bcp)
solver = LinearVariationalSolver(problem)
solver.solve()
(sigma, u, p) = w.split()
# plot
plot(u, axes = True, interactive=True, title = "u")
plot(p, axes = True, interactive=True, title = "p")
Я забыл дх (0) в этом термине. Но это не было проблемой.
В первой части кода (Stokes) я попытался написать условие отсутствия проскальзывания следующим образом:
# Randbedingungen
def top_bottom(x, on_boundary):
return x[1] > 1.0 - DOLFIN_EPS or x[1] < DOLFIN_EPS
noslip = Constant((0.0,0.0))
bc0 = DirichletBC(W.sub(0), noslip, top_bottom)
p_in = 1
p_out = 0
inflow = DirichletBC(W.sub(1), p_in, "x[0] <= 0.0 + DOLFIN_EPS*1000")
outflow = DirichletBC(W.sub(1), p_out, "x[0] >= 0.5 - DOLFIN_EPS*1000")
bcp = [bc0, inflow, outflow]
# Definition f
f = Expression("(0.0, 0.0)")
# Variationsformulierung Stokes
a = inner(grad(u), grad(v))*dx(0) + div(v)*p*dx(0) + q*div(u)*dx(0)
L = inner(f,v)*dx(0)
Но теперь я получаю следующую ошибку:
Shape mismatch: line 56, in <module> L = inner(f,v)*dx(0
Кто-нибудь может помочь? Спасибо!
1 ответ
Я думаю, что здесь есть несколько ошибок. В первой части кода я думаю, что вы используете элементы Тейлора-Гуда для решения уравнения Стокса. Если это так, то U должно быть:
U = VectorFunctionSpace (сетка, "CG", 2)
Также в этой части кода:
a = внутренний (grad (u), grad (v)) * dx + div (v) * p * dx(0) + q * div (u) * dx(0)
L = внутренняя (f,v)*dx(0)
Я не знаю, почему вы не используете dx(0) для первого термина. Я рекомендую вам взглянуть на демонстрации по адресу: http://fenicsproject.org/documentation/dolfin/dev/python/demo/index.html Вы можете получить еще несколько советов.