Coverage for pySDC / projects / StroemungsRaum / problem_classes / NavierStokes_2D_FEniCS.py: 100%

1import logging 

2import dolfin as df 

3import numpy as np 

4import os 

5 

6from pySDC.core.problem import Problem 

7from pySDC.implementations.datatype_classes.fenics_mesh import fenics_mesh, rhs_fenics_mesh 

8 

9 

10class fenics_NSE_2D_mass(Problem): 

11 r""" 

12 Example implementing a forced two-dimensional incompressible Navier–Stokes problem for the DFG 

13 benchmark flow around cylinder using FEniCS (dolfin). 

14 

15 The unknowns are the velocity field :math:`\mathbf{u}(x,y,t)` and the pressure field 

16 :math:`p(x,y,t)` on a 2D domain :math:`\Omega` (here loaded from ``cylinder.xml``). 

17 The model is 

18 

19 .. math:: 

20 \frac{\partial \mathbf{u}}{\partial t} 

21 + (\mathbf{u}\cdot\nabla)\mathbf{u} 

22 - \nu \Delta \mathbf{u} 

23 + \nabla p 

24 = \mathbf{f}, 

25 \qquad 

26 \nabla\cdot \mathbf{u} = 0, 

27 

28 where :math:`\nu` is the kinematic viscosity and :math:`\mathbf{f}` is a forcing term. 

29 

30 This implementation uses a fractional-step /projection method: 

31 a viscous velocity solve, a pressure Poisson solve enforcing incompressibility, and a velocity 

32 correction using the pressure gradient. 

33 

34 Boundary conditions are applied on subsets of the boundary: 

35 - no-slip on channel walls and cylinder surface, 

36 - a time-dependent inflow profile at the inlet, 

37 - pressure condition at the outflow. 

38 

39 Parameters 

40 ---------- 

41 t0 : float, optional 

42 Starting time. 

43 family : str, optional 

44 Finite element family for velocity and pressure spaces. Default is ``'CG'``. 

45 order : int, optional 

46 Polynomial degree for the velocity space. The pressure degree is chosen as ``order - 1``. 

47 nu : float, optional 

48 Kinematic viscosity :math:`\nu`. 

49 

50 Attributes 

51 ---------- 

52 mesh : Mesh 

53 Computational mesh loaded from ``cylinder.xml``. 

54 V : VectorFunctionSpace 

55 Function space for the velocity field. 

56 Q : FunctionSpace 

57 Function space for the pressure field. 

58 M : Matrix 

59 Assembled velocity mass matrix, corresponding to :math:`\int_\Omega \mathbf{u}\cdot\mathbf{v}\,dx`. 

60 K : Matrix 

61 Assembled viscous operator matrix, corresponding to 

62 :math:`- \int_\Omega \nabla \mathbf{u} : (\nu \nabla \mathbf{v})\,dx`. 

63 Mp : Matrix 

64 Assembled pressure mass matrix, corresponding to :math:`\int_\Omega p\,q\,dx`. 

65 g : Expression 

66 Forcing term :math:`\mathbf{f}` (here set to zero). 

67 bcu : list[DirichletBC] 

68 Velocity Dirichlet boundary conditions (walls, inflow, cylinder). 

69 bcp : list[DirichletBC] 

70 Pressure boundary conditions (outflow). 

71 bc_hom : DirichletBC 

72 Homogeneous velocity boundary condition used to fix the residual when requested. 

73 

74 Notes 

75 ----- 

76 The right-hand side evaluation splits into an implicit viscous part and an explicit part 

77 containing forcing and nonlinear convection: 

78 - implicit: :math:`K\,\mathbf{u}` 

79 - explicit: :math:`M\,\mathbf{f} - M\,(\mathbf{u}\cdot\nabla)\mathbf{u}` 

80 """ 

81 

82 dtype_u = fenics_mesh 

83 dtype_f = rhs_fenics_mesh 

84 

85 def __init__(self, t0=0.0, family='CG', order=2, nu=0.001): 

86 

87 # set logger level for FFC and dolfin 

88 logging.getLogger('FFC').setLevel(logging.WARNING) 

89 logging.getLogger('UFL').setLevel(logging.WARNING) 

90 

91 # set solver and form parameters 

92 df.parameters["form_compiler"]["optimize"] = True 

93 df.parameters["form_compiler"]["cpp_optimize"] = True 

94 df.parameters['allow_extrapolation'] = True 

95 

96 # load mesh 

97 path = f"{os.path.dirname(__file__)}/../meshs/" 

98 mesh = df.Mesh(path + '/cylinder.xml') 

99 

100 # define function space for future reference 

101 self.V = df.VectorFunctionSpace(mesh, family, order) 

102 self.Q = df.FunctionSpace(mesh, family, order - 1) 

103 

104 tmp_u = df.Function(self.V) 

105 tmp_p = df.Function(self.Q) 

106 self.logger.debug('Velocity DoFs on this level:', len(tmp_u.vector()[:])) 

107 self.logger.debug('Pressure DoFs on this level:', len(tmp_p.vector()[:])) 

108 

109 super().__init__(self.V) 

110 self._makeAttributeAndRegister('t0', 'family', 'order', 'nu', localVars=locals(), readOnly=True) 

111 

112 # define trial and test functions 

113 self.u = df.TrialFunction(self.V) 

114 self.v = df.TestFunction(self.V) 

115 self.p = df.TrialFunction(self.Q) 

116 self.q = df.TestFunction(self.Q) 

117 

118 # initialize solution functions for pressure 

119 self.pn = df.Function(self.Q) 

120 

121 # mass and stiffness terms 

122 a_K = -1.0 * df.inner(df.nabla_grad(self.u), self.nu * df.nabla_grad(self.v)) * df.dx 

123 a_M = df.inner(self.u, self.v) * df.dx 

124 a_S = df.inner(df.nabla_grad(self.p), df.nabla_grad(self.q)) * df.dx 

125 a_Mp = df.inner(self.p, self.q) * df.dx 

126 

127 # assemble the forms 

128 self.M = df.assemble(a_M) 

129 self.K = df.assemble(a_K) 

130 self.S = df.assemble(a_S) 

131 self.Mp = df.assemble(a_Mp) 

132 

133 # set inflow profile at the domain inlet 

134 self.inflow_profile = df.Expression( 

135 ('4.0*1.5*sin(pi*t/8)*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0'), pi=np.pi, t=t0, degree=self.order 

136 ) 

137 

138 # define boundaries 

139 inflow = 'near(x[0], 0)' 

140 outflow = 'near(x[0], 2.2)' 

141 walls = 'near(x[1], 0) || near(x[1], 0.41)' 

142 cylinder = 'on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3' 

143 

144 # define boundary conditions 

145 bcu_noslip = df.DirichletBC(self.V, df.Constant((0, 0)), walls) 

146 bcu_inflow = df.DirichletBC(self.V, self.inflow_profile, inflow) 

147 bcu_cylinder = df.DirichletBC(self.V, df.Constant((0, 0)), cylinder) 

148 bcp_outflow = df.DirichletBC(self.Q, df.Constant(0), outflow) 

149 

150 # set boundary values 

151 self.bcu = [bcu_noslip, bcu_inflow, bcu_cylinder] 

152 self.bcp = [bcp_outflow] 

153 self.bc_hom = df.DirichletBC(self.V, df.Constant((0.0, 0.0)), 'on_boundary') 

154 self.fix_bc_for_residual = True 

155 

156 # Define measure for drag and lift computation 

157 Cylinder = df.CompiledSubDomain(cylinder) 

158 CylinderBoundary = df.MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0) 

159 Cylinder.mark(CylinderBoundary, 1) 

160 self.dsc = df.Measure("ds", domain=mesh, subdomain_data=CylinderBoundary, subdomain_id=1) 

161 

162 # set forcing term as expression 

163 self.g = df.Expression(('0.0', '0.0'), t=self.t0, degree=self.order) 

164 

165 # initialize XDMF files for velocity and pressure if needed 

166 self.xdmffile_p = None 

167 self.xdmffile_u = None 

168 

169 def solve_system(self, rhs, factor, u0, t, dtau): 

170 r""" 

171 Dolfin's linear solver for :math:`(M - factor A) \vec{u} = \vec{rhs}`. 

172 

173 Parameters 

174 ---------- 

175 rhs : dtype_f 

176 Right-hand side for the nonlinear system. 

177 factor : float 

178 Factor for the linear system. 

179 u0 : dtype_u 

180 Initial guess for the iterative solver (not used here so far). 

181 t : float 

182 Current time. 

183 dtau : float 

184 Zero-to-node step size in the SDC sweep, used in the pressure correction step. 

185 

186 Returns 

187 ------- 

188 u : dtype_u 

189 Solution: velocity. 

190 p : space function 

191 Solution: pressure. 

192 

193 """ 

194 

195 u = self.dtype_u(u0) 

196 p = df.Function(self.Q) 

197 T = self.M - factor * self.K 

198 

199 # solve for the intermediate velocity 

200 self.inflow_profile.t = t 

201 [bc.apply(T, rhs.values.vector()) for bc in self.bcu] 

202 df.solve(T, u.values.vector(), rhs.values.vector()) 

203 

204 # solve for the pressure correction 

205 L2 = -(1 / dtau) * df.div(u.values) * self.q * df.dx 

206 Rhs2 = df.assemble(L2) 

207 [bc.apply(self.S, Rhs2) for bc in self.bcp] 

208 df.solve(self.S, p.vector(), Rhs2) 

209 

210 # velocity correction 

211 L3 = df.dot(u.values, self.v) * df.dx - dtau * df.dot(df.nabla_grad(p), self.v) * df.dx 

212 Rhs3 = df.assemble(L3) 

213 [bc.apply(self.M, Rhs3) for bc in self.bcu] 

214 df.solve(self.M, u.values.vector(), Rhs3) 

215 

216 return u, p 

217 

218 def eval_f(self, u, t): 

219 """ 

220 Routine to evaluate both parts of the right-hand side. 

221 

222 Parameters 

223 ---------- 

224 u : dtype_u 

225 Current values of the numerical solution. 

226 t : float 

227 Current time at which the numerical solution is computed. 

228 

229 Returns 

230 ------- 

231 f : dtype_f 

232 The right-hand side divided into two parts. 

233 """ 

234 self.g.t = t 

235 conv = -1.0 * df.dot(u.values, df.nabla_grad(u.values)) 

236 

237 f = self.dtype_f(self.V) 

238 

239 self.K.mult(u.values.vector(), f.impl.values.vector()) 

240 f.expl = self.apply_mass_matrix(self.dtype_u(df.project(conv, self.V))) 

241 f.expl += self.apply_mass_matrix(self.dtype_u(df.interpolate(self.g, self.V))) 

242 

243 return f 

244 

245 def apply_mass_matrix(self, u): 

246 r""" 

247 Routine to apply mass matrix. 

248 

249 Parameters 

250 ---------- 

251 u : dtype_u 

252 Current values of the numerical solution. 

253 

254 Returns 

255 ------- 

256 me : dtype_u 

257 The product :math:`M \vec{u}`. 

258 """ 

259 

260 me = self.dtype_u(self.V) 

261 self.M.mult(u.values.vector(), me.values.vector()) 

262 

263 return me 

264 

265 def u_exact(self, t): 

266 r""" 

267 Routine to compute the exact solution at time :math:`t`. 

268 

269 Parameters 

270 ---------- 

271 t : float 

272 Time of the exact solution. 

273 

274 Returns 

275 ------- 

276 me : dtype_u 

277 Exact solution. 

278 """ 

279 

280 u0 = df.Expression(('0.0', '0.0'), degree=self.order) 

281 me = self.dtype_u(df.interpolate(u0, self.V), val=self.V) 

282 

283 return me 

284 

285 def fix_residual(self, res): 

286 """ 

287 Applies homogeneous Dirichlet boundary conditions to the residual 

288 

289 Parameters 

290 ---------- 

291 res : dtype_u 

292 Residual 

293 """ 

294 self.bc_hom.apply(res.values.vector()) 

295 return None