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

1import logging 

2import os 

3 

4import dolfin as df 

5import numpy as np 

6 

7from pySDC.core.problem import Problem 

8from pySDC.implementations.datatype_classes.fenics_mesh import fenics_mesh 

9from pySDC.projects.StroemungsRaum.problem_classes.newton_step import NewtonStep 

10 

11 

12class fenics_NSE_2D_Monolithic(Problem): 

13 r""" 

14 Example implementing the forced two-dimensional incompressible Navier-Stokes equations with 

15 time-dependent Dirichlet boundary conditions 

16 

17 .. math:: 

18 \frac{\partial u}{\partial t} = - u \cdot \nabla u + \nu \Delta u - \nabla p + f 

19 0 = \nabla \cdot u 

20 

21 for :math:`x \in \Omega`, where the forcing term :math:`f` is defined by 

22 

23 .. math:: 

24 f(x, t) = (0, 0). 

25 

26 This implementation follows a monolithic approach, where velocity and pressure are 

27 solved simultaneously in a coupled system using a mixed finite element formulation. 

28 

29 Boundary conditions are applied on subsets of the boundary: 

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

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

32 - pressure condition at the outflow. 

33 

34 In this class the problem is implemented in the way that the spatial part is solved using ``FEniCS`` [1]_. Hence, the problem 

35 is reformulated to the *weak formulation* 

36 

37 .. math:: 

38 \int_\Omega u_t v\,dx = - \int_\Omega u \cdot \nabla u v\,dx - \nu \int_\Omega \nabla u \nabla v\,dx + \int_\Omega p \nabla \cdot v\,dx + \int_\Omega f v\,dx 

39 \int_\Omega \nabla \cdot u q\,dx = 0 

40 

41 Parameters 

42 ---------- 

43 t0 : float, optional 

44 Starting time. 

45 order : int, optional 

46 Defines the order of the elements in the function space. 

47 nu : float, optional 

48 Diffusion coefficient :math:`\nu`. 

49 Sol_tol : float, optional 

50 Tolerance for the nonlinear solver. 

51 

52 Attributes 

53 ---------- 

54 V : FunctionSpace 

55 Defines the function space of the trial and test functions for velocity. 

56 Q : FunctionSpace 

57 Defines the function space of the trial and test functions for pressure. 

58 W : FunctionSpace 

59 Defines the mixed function space for the coupled velocity-pressure system. 

60 M : scalar, vector, matrix or higher rank tensor 

61 Denotes the expression :math:`\int_\Omega u_t v\,dx`. 

62 g : Expression 

63 The forcing term :math:`f` in the Navier-Stokes momentum equation. 

64 bc : DirichletBC 

65 Denotes the time-dependent Dirichlet boundary conditions. 

66 bc_hom : DirichletBC 

67 Denotes the homogeneous Dirichlet boundary conditions, potentially required for fixing the residual 

68 fix_bc_for_residual: boolean 

69 flag to indicate that the residual requires special treatment due to boundary conditions 

70 

71 References 

72 ---------- 

73 .. [1] The FEniCS Project Version 1.5. M. S. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, 

74 C. Richardson, J. Ring, M. E. Rognes, G. N. Wells. Archive of Numerical Software (2015). 

75 .. [2] Automated Solution of Differential Equations by the Finite Element Method. A. Logg, K.-A. Mardal, G. N. 

76 Wells and others. Springer (2012). 

77 """ 

78 

79 dtype_u = fenics_mesh 

80 dtype_f = fenics_mesh 

81 

82 df.set_log_active(False) 

83 

84 def __init__(self, t0=0.0, order=2, nu=0.001, Sol_tol=1e-10): 

85 

86 # set logger level for FFC and dolfin 

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

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

89 

90 # set solver and form parameters 

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

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

93 df.parameters['allow_extrapolation'] = True 

94 

95 # load mesh 

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

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

98 

99 # define function spaces for future reference (Taylor-Hood) 

100 P2 = df.VectorElement("P", mesh.ufl_cell(), order) 

101 P1 = df.FiniteElement("P", mesh.ufl_cell(), order - 1) 

102 TH = df.MixedElement([P2, P1]) 

103 self.W = df.FunctionSpace(mesh, TH) 

104 self.V = df.FunctionSpace(mesh, P2) 

105 self.Q = df.FunctionSpace(mesh, P1) 

106 

107 # print the number of DoFs for debugging purposes 

108 tmp = df.Function(self.W) 

109 self.logger.debug('DoFs on this level: %d', len(tmp.vector()[:])) 

110 

111 super().__init__(self.W) 

112 self._makeAttributeAndRegister('t0', 'order', 'nu', 'Sol_tol', localVars=locals(), readOnly=True) 

113 

114 # Trial and test function for the Mixed FE space 

115 self.u, self.p = df.TrialFunctions(self.W) 

116 self.v, self.q = df.TestFunctions(self.W) 

117 

118 # velocity mass matrix 

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

120 self.M = df.assemble(a_M) 

121 

122 # define the time-dependent inflow profile as an Expression 

123 Uin = '4.0*1.5*sin(pi*t/8)*x[1]*(0.41 - x[1]) / pow(0.41, 2)' 

124 self.u_in = df.Expression((Uin, '0'), pi=np.pi, t=t0, degree=self.order) 

125 

126 # define boundaries 

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

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

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

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

131 

132 # define boundary conditions 

133 bc_in = df.DirichletBC(self.W.sub(0), self.u_in, inflow) 

134 bc_out = df.DirichletBC(self.W.sub(1), 0, outflow) 

135 bc_walls = df.DirichletBC(self.W.sub(0), (0, 0), walls) 

136 bc_cylinder = df.DirichletBC(self.W.sub(0), (0, 0), cylinder) 

137 self.bc = [bc_cylinder, bc_walls, bc_out, bc_in] 

138 

139 # homogeneous boundary conditions for fixing the residual 

140 bc_hom_u = df.DirichletBC(self.W.sub(0), df.Constant((0, 0)), 'on_boundary') 

141 bc_hom_p = df.DirichletBC(self.W.sub(1), df.Constant(0), 'on_boundary') 

142 self.bc_hom = [bc_hom_u, bc_hom_p] 

143 self.fix_bc_for_residual = True 

144 

145 # define measure for drag and lift computation 

146 Cylinder = df.CompiledSubDomain(cylinder) 

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

148 Cylinder.mark(CylinderBoundary, 1) 

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

150 

151 # set forcing term as expression 

152 self.g = df.Expression(('0', '0'), a=np.pi, b=self.nu, t=self.t0, degree=self.order) 

153 

154 # initialize XDMF files for velocity and pressure if needed 

155 self.xdmffile_p = None 

156 self.xdmffile_u = None 

157 

158 # interpolated forcing term referenced by the residual form 

159 self.g_h = df.Function(self.V) 

160 

161 # residual form for a single node-to-node step, assembled once; the step size is a 

162 # Constant and the boundary and forcing data carry the time dependence 

163 self.factor = df.Constant(0.0) 

164 self.w = df.Function(self.W) 

165 u, p = df.split(self.w) 

166 

167 F = df.dot(u, self.v) * df.dx 

168 F += self.factor * df.dot(df.dot(u, df.nabla_grad(u)), self.v) * df.dx 

169 F += self.factor * self.nu * df.inner(df.nabla_grad(u), df.nabla_grad(self.v)) * df.dx 

170 F -= self.factor * df.dot(p, df.div(self.v)) * df.dx 

171 F -= self.factor * df.dot(self.g_h, self.v) * df.dx 

172 F -= self.factor * df.dot(df.div(u), self.q) * df.dx 

173 

174 self.step = NewtonStep(F, df.derivative(F, self.w)) 

175 self.newton = df.NewtonSolver() 

176 self.newton.parameters['absolute_tolerance'] = Sol_tol 

177 

178 def use_restricted_boundary_values(self, u, t): 

179 """ 

180 Refuse multilevel runs: the solves impose the time-dependent boundary data, which at the coarse nodes is not 

181 the restricted fine boundary values, so MLSDC and PFASST would not converge to the collocation solution. See 

182 :meth:`pySDC.core.problem.Problem.use_restricted_boundary_values`. 

183 """ 

184 raise NotImplementedError(f'{type(self).__name__} cannot impose restricted boundary values on coarse levels') 

185 

186 def solve_system(self, rhs, factor, u0, t): 

187 r""" 

188 Dolfin's nonlinear solver for : 

189 .. math:: 

190 (u,v) + factor * (u \cdot \nabla u, v) + factor * \nu (\nabla u, \nabla v) - factor * (p, \nabla \cdot v) - factor * (g, v) - factor * (div(u), q) = (rhs_u, v) + (rhs_p, q) 

191 

192 Parameters 

193 ---------- 

194 rhs : dtype_f 

195 Right-hand side for the nonlinear system. 

196 factor : float 

197 Abbrev. for the node-to-node stepsize (or any other factor required). 

198 u0 : dtype_u 

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

200 t : float 

201 Current time. 

202 

203 Returns 

204 ------- 

205 w : dtype_u 

206 Solution. 

207 """ 

208 # update time in boundary conditions and in the forcing term 

209 self.u_in.t = t 

210 self.g.t = t 

211 self.g_h.interpolate(self.g) 

212 self.factor.assign(factor) 

213 

214 # the SDC right-hand side enters the residual as a vector, no mass matrix involved 

215 self.w.vector()[:] = u0.values.vector()[:] 

216 self.step.rhs = rhs.values.vector() 

217 self.step.bcs = self.bc 

218 

219 # solve the nonlinear system using Newton's method 

220 self.newton.solve(self.step, self.w.vector()) 

221 

222 w = self.dtype_u(self.W) 

223 w.values.vector()[:] = self.w.vector()[:] 

224 

225 return w 

226 

227 def eval_f(self, w, t): 

228 """ 

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

230 

231 Parameters 

232 ---------- 

233 w : dtype_u 

234 Current values of the numerical solution. 

235 t : float 

236 Current time at which the numerical solution is computed. 

237 

238 Returns 

239 ------- 

240 f : dtype_f 

241 The right-hand side. 

242 """ 

243 

244 f = self.dtype_f(self.W) 

245 

246 u, p = df.split(w.values) 

247 

248 # get the forcing term 

249 self.g.t = t 

250 g = self.dtype_f(df.interpolate(self.g, self.V)) 

251 

252 F = -1.0 * df.dot(df.dot(u, df.nabla_grad(u)), self.v) * df.dx 

253 F -= self.nu * df.inner(df.nabla_grad(u), df.nabla_grad(self.v)) * df.dx 

254 F += df.dot(p, df.div(self.v)) * df.dx 

255 F += df.dot(g.values, self.v) * df.dx 

256 F += df.dot(df.div(u), self.q) * df.dx 

257 

258 df.assemble(F, tensor=f.values.vector()) 

259 

260 return f 

261 

262 def apply_mass_matrix(self, w): 

263 r""" 

264 Routine to apply velocity mass matrix. 

265 

266 Parameters 

267 ---------- 

268 w : dtype_u 

269 Current values of the numerical solution. 

270 

271 Returns 

272 ------- 

273 me : dtype_u 

274 The product :math:`M \vec{w}`. 

275 """ 

276 

277 me = self.dtype_u(self.W) 

278 self.M.mult(w.values.vector(), me.values.vector()) 

279 

280 return me 

281 

282 def u_exact(self, t): 

283 r""" 

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

285 

286 Parameters 

287 ---------- 

288 t : float 

289 Time of the exact solution. 

290 

291 Returns 

292 ------- 

293 me : dtype_u 

294 Exact solution. 

295 """ 

296 # define the exact solution 

297 w = df.Function(self.W) 

298 

299 # assign the exact solution as an Expression 

300 df.assign(w.sub(0), df.interpolate(df.Expression(('0.0', '0.0'), degree=self.order), self.V)) 

301 df.assign(w.sub(1), df.interpolate(df.Expression('0.0', degree=self.order - 1), self.Q)) 

302 

303 # update time in Boundary conditions 

304 self.u_in.t = t 

305 [bc.apply(w.vector()) for bc in self.bc] 

306 

307 me = self.dtype_u(w, val=self.W) 

308 

309 return me 

310 

311 def fix_residual(self, res): 

312 """ 

313 Applies homogeneous Dirichlet boundary conditions to the residual 

314 

315 Parameters 

316 ---------- 

317 res : dtype_u 

318 Residual 

319 """ 

320 [bc.apply(res.values.vector()) for bc in self.bc_hom] 

321 

322 return None