firedrake创建mesh和求解具体问题

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

界面问题

$-\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\ge 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}$
$u=h,x\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 vdx$
$b={\int }_{{\Gamma }_{in}}{g}_i\ast vds+{\int }_{{\Gamma }_{out}}{g}_{out}\ast vds+{\int }_{\Omega }f\ast v\ast dx$
在firedrake中需要设置把对应的区域和不同的$\Gamma$表示出来，有了上面的ID以后就有
$a={\int }_{\Omega }\nabla u\cdot \nabla vdx$
$b={\int }_{{\Gamma }_{in}}{g}_i\ast vds(17)+{\int }_{{\Gamma }_{out}}{g}_{out}\ast vds(18)+{\int }_{\Omega }f\ast v\ast 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).
$\{\begin{array}{l}-\nu \Delta u+u{u}_x+v{u}_y+p_x=f_1,\\ -\nu \Delta v+u{v}_x+v{v}_y+p_y=f_2,\\ u_x+v_y=0,\\ u{|}_{\partial \Omega }=u_0,\\ v{|}_{\partial \Omega }=v_0.\end{array}$
$\begin{array}{rl}{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{array}$
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 --