firedrake创建mesh和求解具体问题

Galerkin码农选手

已于 2023-04-12 19:40:52 修改

阅读量435

点赞数

分类专栏：传统数值方法 NS方程求解文章标签： mesh

于 2022-05-06 10:39:58 首次发布

本文链接：https://blog.csdn.net/forrestguang/article/details/124552846

版权

传统数值方法同时被 2 个专栏收录

35 篇文章 20 订阅 ¥199.90 ¥299.90

订阅专栏

NS方程求解

8 篇文章 28 订阅

订阅专栏

这篇博客介绍了如何使用firedrake解决界面问题和bump泊松问题。首先，详细解释了如何通过Gmesh创建mesh，并设定边界条件。接着，展示了在bump.geo文件中定义不同区域的过程。然后，讨论了泊松方程的双线性形式，并解释了在firedrake中如何处理Dirichlet边界条件。最后，提到了kflow问题的求解，以及如何通过paraview查看和分析解的结果。

摘要由CSDN通过智能技术生成

博客内容参考知乎添加链接描述

界面问题

$-\nabla \cdot(\sigma \nabla u) + u = 5,\Omega$
$u|_{\partial \Omega_1} = 0$
$[\sigma\nabla u \cdot n] = 3,in \partial \Omega_2$
$\partial \Omega_1$ 为矩形边界
$\partial \Omega_2$ 为圆形边界
在这里插入图片描述
geo文件和代码参考链接添加链接描述
使用paraview得到的结果如下，运行下面的interface.py，得到face_0.vtu文件，利用paraview打开该文件即可得到

test1.geo

Point(1) = {-6,  2, 0, 0.5};
Point(2) = {-6, -2, 0, 0.5};
Point(3) = { 6, -2, 0, 0.5};
Point(4) = { 6,  2, 0, 0.5};
Point(5) = { 0,  0, 0, 0.1};
Point(6) = { 1,  0, 0, 0.1};
Point(7) = {-1,  0, 0, 0.1};
Point(8) = { 0,  1, 0, 0.1};
Point(9) = { 0, -1, 0, 0.1};
Line(1) = {1, 4};
Line(2) = {4, 3};
Line(3) = {3, 2};
Line(4) = {2, 1};
Circle(5) = {8, 5, 6};
Circle(6) = {6, 5, 9};
Circle(7) = {9, 5, 7};
Circle(8) = {7, 5, 8};
Curve Loop( 9) = {1, 2, 3, 4};
Curve Loop(10) = {8, 5, 6, 7};
Plane Surface(1) = {9, 10};
Plane Surface(2) = {10};
Physical Curve("HorEdges", 11) = {1, 3};
Physical Curve("VerEdges", 12) = {2, 4};
Physical Curve("Circle", 13) = {8,5,6,7};
Physical Surface("PunchedDom", 3) = {1};
Physical Surface("Disc", 4) = {2};

interface.py

from firedrake import *

sigma1 = 1
sigma2 = 2
#mesh = Mesh('immersed_domain.msh')
mesh = Mesh('test1.msh')
V = FunctionSpace(mesh,'CG',1)

u = TrialFunction(V)
v = TestFunction(V)
a = sigma2*dot(grad(v),grad(u))*dx(4) + sigma1*dot(grad(u),grad(v))*dx(3) + u*v*dx
L = Constant(5.)*v*dx + Constant(3.)*v('+')*dS(13)

Dirbc = DirichletBC(V,0,[11,12])

u = Function(V)
solve(a == L,u,bcs = Dirbc,solver_parameters = {'ksp_type':'cg'})
File('face.pvd').write(u)

bump泊松问题

现在考虑区域
在这里插入图片描述
其中左右分别是inflow,outflow,上面的圆弧是free,其余都是wall，
$\Omega=[-1,1]\times[-1,1]\cup\Omega_1$
$\Omega_1=\{(x,y):x^2 + (y-1)^2 < 0.5^2,y \geq 1\}$
首先需要用Gmesh生成msh文件，代码如下：
对于其中的Point,Line,Circle和Curve Loop，读者只需要查看最简单的Gmesh入门即可，这里需要注意的是最后几行Physical开头的命令，如果这几行命令不写，firedrake无法识别生成的msh文件，缺乏边界条件，这几行是用来告诉firedrake边界条件的，最后一行Physical Surface如果缺失，firedrake无法读取msh文件，如果只写Physical Curve这几行，firedrake读取的msh文件会发现得到的是一维剖分
其中Physical Curve(‘xx’,10)={1,2};表示把Line(1),Line(2)标记为xx，ID是10 .具体参考本人的geo文件。

bump.geo

lr1 = 5e-1;
lr2 = 1e-1;

Point(1) = {-1,1,0,lr1};
Point(2) = {-1,-1,0,lr1};
Point(3) = {1,-1,0,lr1};
Point(4) = {1,1,0,lr1};

Point(5) = {0.5,1,0,lr2};
Point(6) = {0,1,0,lr2};
Point(7) = {0,1.5,0,lr2};
Point(8) = {-0.5,1,0,lr2};

Line(9) = {1,2};
Line(10) = {2,3};
Line(11) = {3,4};
Line(12) = {4,5};
Circle(13) = {5,6,7};
Circle(14) = {7,6,8};

Line(15) = {8,1};

Curve Loop(16) = {9,10,11,12,13,14,15};

Plane Surface(1) = {16};

Physical Curve("inflow",17) = {9};
Physical Curve("outflow",18) = {11};
Physical Curve("wall",19) = {10,12,15};
Physical Curve("free",20) = {13,14};

Physical Surface("PunchedDom",21) = {1};

先求解一个简单的泊松方程来做实验：
考虑问题
$-\Delta u=f,x \in \Omega$
$\frac{\partial u}{\partial n}=g_i,x \in \Gamma_{in}$
$\frac{\partial u}{\partial n}=g_o,x \in \Gamma_{out}$
$\in \Gamma_{wall}\cup\Gamma_{free}$
我们在上面的bump.geo文件中设定
Physical Curve(“inflow”,17) = {9};这一行表示把Line(9)定义inflow,用ID=17表示，其他类似
Physical Curve(“outflow”,18) = {11};
Physical Curve(“wall”,19) = {10,12,15};这一行表示把Line(10)，Line（12），Line（15）定义wall,用ID=19表示，
Physical Curve(“free”,20) = {13,14};
双线性形式如下：
$a=\int_{\Omega}\nabla u \cdot \nabla v dx$
$\int_{\Gamma_{in}}g_i*vds + \int_{\Gamma_{out}}g_{out}*vds + \int_{\Omega} f*v*dx$
在firedrake中需要设置把对应的区域和不同的 $\Gamma$ 表示出来，有了上面的ID以后就有
$a=\int_{\Omega}\nabla u \cdot \nabla v dx$
$\int_{\Gamma_{in}}g_i*vds(17) + \int_{\Gamma_{out}}g_{out}*vds(18) + \int_{\Omega} f*v*dx$
我们还需要表示剩下的Dirchlet边界条件，有了上面的ID，我们就可以得到
bc = DirichletBC(V,h,[19,20])，代入到求解函数中即可。

bumkpoi.py

在这里插入图片描述

from firedrake import *

mesh = Mesh('bump.msh')
x = SpatialCoordinate(mesh)


V = FunctionSpace(mesh,'CG',1)
#print(mesh.coordinates.dat.data)
u = TrialFunction(V)
v = TestFunction(V)
f = Function(V)

f.interpolate((pi*pi - 1)*sin(pi*x[0])*exp(x[1]))
g_o = Function(V).interpolate(-pi*exp(x[1]))
g_i = Function(V).interpolate(pi*exp(x[1]))
h = Function(V)
h.interpolate(exp(x[1])*sin(pi*x[0]))
a = dot(grad(u),grad(v))*dx

L = f*v*dx + g_i*v*ds(17) + g_o*v*ds(18)
DirBC = DirichletBC(V,h,[19,20])
u = Function(V)
solve(a == L, u, bcs=DirBC, solver_parameters={'ksp_type': 'cg'})
u_acc = Function(V).interpolate(exp(x[1])*sin(pi*x[0]))
L1_err = assemble(abs(u - u_acc)*dx)
print('firedrake solve poisson:L1_err:%3e'%(L1_err))
File('bump.pvd').write(u)

bumpkflow.py

考虑上述区域的kflow问题
参考文献：Dockhorn, T. . “A Discussion on Solving Partial Differential Equations using Neural Networks.” (2019).
$\left\{\begin{array}{l} -\nu \Delta u + uu_x + v u_y + p_x = f_1,\\ -\nu \Delta v + uv_x + vv_y + p_y=f_2,\\ u_x + v_y=0,\\ u|_{\partial \Omega} = u_0,\\ v|_{\partial \Omega}=v_0. \end{array}\right.$
$\begin{aligned} u_{1}&\left(x_{1}, x_{2}\right) =1-\exp \left(\lambda x_{1}\right) \cos 2 \pi x_{2}, \\ u_{2}&\left(x_{1}, x_{2}\right) =\frac{\lambda}{2 \pi} \exp \left(\lambda x_{1}\right) \sin 2 \pi x_{2}, \\ p&\left(x_{1}, x_{2}\right) =\frac{1}{2}\left(1-\exp \left(2 \lambda x_{1}\right)\right)+C,\\ \lambda&=\frac{1}{2\nu} - \sqrt{\frac{1}{4\nu^2} + 4\pi^2},\nu = 0.025,\\ f_1&=0,f_2=0. \end{aligned}$
paraview使用，我们在代码里面写了
u,p = solution.split()
u.rename(“kflowu”)
p.rename(“kflowp”)
File(“kflow.pvd”).write(u,p)
这四行代码，会生成kflow.pvd和kflow_0.vtu，paraview加载kflow_0.vtu文件，一开始可能是solid color，修改红色圈里面的参数pred[0]就可以正常显示解

在这里插入图片描述

可以通过refine mesh来查看加细之前和加细之后的误差比例，如果误差满足一个order的减少，说明代码就没问题

from firedrake import *
import numpy as np

Re = 40
nu = 1/Re
mesh = Mesh('bump.msh')
x = SpatialCoordinate(mesh)


V = VectorFunctionSpace(mesh,'CG',2)*FunctionSpace(mesh,'CG',1)
solution = Function(V,name = 'pred')
testfunction = TestFunction(V)
u,p = split(solution)
v,q = split(testfunction)

lam = Constant(1/(2*nu) - sqrt(1/(4*nu*nu) + 4*pi*pi))
print(1/(2*nu) - sqrt(1/(4*nu*nu) + 4*pi*pi))
u_acc = as_vector([1 - exp(lam*x[0])*cos(2*pi*x[1]),lam*exp(lam*x[0])*sin(2*pi*x[1])/(2*pi)])
#print(mesh.coordinates.dat.data)

a = nu*inner(grad(u),grad(v))*dx + inner(dot(grad(u),u),v)*dx - p*div(v)*dx + div(u)*q*dx

DirBC = DirichletBC(V.sub(0),u_acc,[17,18,19,20])
parameter = {"snes_monitor": None,
                             "snes_rtol":1e-4,"snes_atol":1e-50,"snes_stol":1e-8,"snes_max_it":100,
                             "ksp_type": "gmres",
                             "mat_type": "aij",
                             "pc_type": "lu",
                             "pc_factor_mat_solver_type": "mumps"}
#parameter = {'ksp_type':'cg'}
solve(a == 0, solution, bcs=DirBC, solver_parameters=parameter)



x = mesh.coordinates.dat.data
u,p = solution.split()
u.rename("kflowu")
p.rename("kflowp")
File("kflow.pvd").write(u,p)
import matplotlib.pyplot as plt
fig, axes = plt.subplots()
levels = np.linspace(0, 1, 51)
contours = tricontourf(p, levels=levels, axes=axes, cmap="inferno")
axes.set_aspect("equal")
fig.colorbar(contours)
fig.show()
plt.savefig('kflow.png')
err = inner(u - u_acc,u - u_acc)*dx
print(err)
L1_err = assemble(abs(u[0] - u_acc[0])*dx)
print('the obj:%.3e'%(L1_err))
-- INSERT --