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

1import logging 

2 

3import dolfin as df 

4import numpy as np 

5 

6from pySDC.core.problem import Problem 

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

8 

9 

10class fenics_heat2D_mass(Problem): 

11 r""" 

12 Example implementing the forced two-dimensional heat equation with Dirichlet boundary conditions 

13 

14 .. math:: 

15 \frac{\partial u}{\partial t} = \nu \Delta u + f 

16 

17 for :math:`x \in \Omega:=[0,1]x[0,1]`, where the forcing term :math:`f` is defined by 

18 

19 .. math:: 

20 f(x, y, t) = -\sin(\pi x)\sin(\pi y) (\sin(t) - 2\nu \pi^2 \cos(t)). 

21 

22 For initial conditions with constant c and 

23 

24 .. math:: 

25 u(x, y, 0) = \sin(\pi x)\sin(\pi y) + c 

26 

27 the exact solution of the problem is given by 

28 

29 .. math:: 

30 u(x, y, t) = \sin(\pi x)\sin(\pi y)\cos(t) + c. 

31 

32 In this class the spatial part is solved using ``FEniCS`` [1]_. Hence, the problem is reformulated to 

33 the *weak formulation* 

34 

35 .. math: 

36 \int_\Omega u_t v\,dx = - \nu \int_\Omega \nabla u \nabla v\,dx + \int_\Omega f v\,dx. 

37 

38 The part containing the forcing term is treated explicitly, where it is interpolated in the function 

39 space. The other part will be treated in an implicit way. 

40 

41 Parameters 

42 ---------- 

43 c_nvars : int, optional 

44 Spatial resolution. 

45 t0 : float, optional 

46 Starting time. 

47 family : str, optional 

48 Indicates the family of elements used to create the function space 

49 for the trail and test functions. The default is ``'CG'``, which are the class 

50 of Continuous Galerkin, a *synonym* for the Lagrange family of elements, see [2]_. 

51 order : int, optional 

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

53 nu : float, optional 

54 Diffusion coefficient :math:`\nu`. 

55 c: float, optional 

56 Constant for the Dirichlet boundary condition :math: `c` 

57 

58 Attributes 

59 ---------- 

60 V : FunctionSpace 

61 Defines the function space of the trial and test functions. 

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

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

64 K : scalar, vector, matrix or higher rank tensor 

65 Denotes the expression :math:`- \nu \int_\Omega \nabla u \nabla v\,dx`. 

66 g : Expression 

67 The forcing term :math:`f` in the heat equation. 

68 bc : DirichletBC 

69 Denotes the Dirichlet 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 = rhs_fenics_mesh 

81 

82 def __init__(self, c_nvars=64, t0=0.0, family='CG', order=2, nu=0.1, c=0.0): 

83 

84 # define the Dirichlet boundary 

85 def Boundary(x, on_boundary): 

86 return on_boundary 

87 

88 # set solver and form parameters 

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

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

91 df.parameters['allow_extrapolation'] = True 

92 

93 # set mesh and refinement (for multilevel) 

94 mesh = df.UnitSquareMesh(c_nvars, c_nvars) 

95 

96 # define function space for future reference 

97 self.V = df.FunctionSpace(mesh, family, order) 

98 tmp = df.Function(self.V) 

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

100 

101 # invoke super init, passing number of dofs, dtype_u and dtype_f 

102 super().__init__(self.V) 

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

104 

105 # Stiffness term (Laplace) 

106 u = df.TrialFunction(self.V) 

107 v = df.TestFunction(self.V) 

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

109 

110 # Mass term 

111 a_M = u * v * df.dx 

112 

113 self.M = df.assemble(a_M) 

114 self.K = df.assemble(a_K) 

115 

116 # set boundary values 

117 self.u_D = df.Expression('sin(a*x[0]) * sin(a*x[1]) * cos(t) + c', c=self.c, a=np.pi, t=t0, degree=self.order) 

118 self.bc = df.DirichletBC(self.V, self.u_D, Boundary) 

119 self.bc_hom = df.DirichletBC(self.V, df.Constant(0), Boundary) 

120 

121 self.fix_bc_for_residual = True 

122 

123 # set forcing term as expression 

124 self.g = df.Expression( 

125 '-sin(a*x[0]) * sin(a*x[1]) * (sin(t) - 2*b*a*a*cos(t))', 

126 a=np.pi, 

127 b=self.nu, 

128 t=self.t0, 

129 degree=self.order, 

130 ) 

131 

132 def use_restricted_boundary_values(self, u, t): 

133 """ 

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

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

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

137 """ 

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

139 

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

141 r""" 

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

143 

144 Parameters 

145 ---------- 

146 rhs : dtype_f 

147 Right-hand side for the nonlinear system. 

148 factor : float 

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

150 u0 : dtype_u 

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

152 t : float 

153 Current time. 

154 

155 Returns 

156 ------- 

157 u : dtype_u 

158 Solution. 

159 """ 

160 

161 u = self.dtype_u(u0) 

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

163 b = self.dtype_u(rhs) 

164 

165 self.u_D.t = t 

166 

167 self.bc.apply(T, b.values.vector()) 

168 self.bc.apply(b.values.vector()) 

169 

170 df.solve(T, u.values.vector(), b.values.vector()) 

171 

172 return u 

173 

174 def eval_f(self, u, t): 

175 """ 

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

177 

178 Parameters 

179 ---------- 

180 u : dtype_u 

181 Current values of the numerical solution. 

182 t : float 

183 Current time at which the numerical solution is computed. 

184 

185 Returns 

186 ------- 

187 f : dtype_f 

188 The right-hand side divided into two parts. 

189 """ 

190 

191 f = self.dtype_f(self.V) 

192 

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

194 

195 self.g.t = t 

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

197 return f 

198 

199 def apply_mass_matrix(self, u): 

200 r""" 

201 Routine to apply mass matrix. 

202 

203 Parameters 

204 ---------- 

205 u : dtype_u 

206 Current values of the numerical solution. 

207 

208 Returns 

209 ------- 

210 me : dtype_u 

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

212 """ 

213 

214 me = self.dtype_u(self.V) 

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

216 

217 return me 

218 

219 def u_exact(self, t): 

220 r""" 

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

222 

223 Parameters 

224 ---------- 

225 t : float 

226 Time of the exact solution. 

227 

228 Returns 

229 ------- 

230 me : dtype_u 

231 Exact solution. 

232 """ 

233 u0 = df.Expression('sin(a*x[0]) * sin(a*x[1]) * cos(t) + c', c=self.c, a=np.pi, t=t, degree=self.order) 

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

235 

236 return me 

237 

238 def fix_residual(self, res): 

239 """ 

240 Applies homogeneous Dirichlet boundary conditions to the residual 

241 

242 Parameters 

243 ---------- 

244 res : dtype_u 

245 Residual 

246 """ 

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

248 return None