FEniCS:UMFPACK报告说正在求解的矩阵是singu

2024-05-15 08:05:47 发布

您现在位置:Python中文网/ 问答频道 /正文

我创建了一个分区域,第一个子域中的stokes方程和第二个子域中的混合poisson方程(darcy)。我使用UnitSquare,子域1应该是0到0,5的区间,子域2应该是0,5到1的区间。在

但现在我得到了以下错误:

求解线性变分问题。 与调用numeric相关的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")

我忘了这学期的dx(0)。但这不是问题所在。在

在代码的第一部分(Stokes)中,我尝试用以下方式编写防滑条件:

^{pr2}$

但现在我得到了以下错误:

Shape mismatch: line 56, in <module> L = inner(f,v)*dx(0

有人能帮忙吗?谢谢!在


Tags: 子域inepscgdesmeshdolfinboundary
1条回答
网友
1楼 · 发布于 2024-05-15 08:05:47

我认为这里有几个错误。在代码的第一部分,我想你是用泰勒胡德元素来解斯托克斯方程。如果是这样的话,那么你应该是:

U=向量函数空间(网格,“CG”,2)

在本部分代码中:

a=内部(梯度(u),梯度(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 你可能会得到更多的提示。在

相关问题 更多 >