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 solve_system(self, rhs, factor, u0, t): 

133 r""" 

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

135 

136 Parameters 

137 ---------- 

138 rhs : dtype_f 

139 Right-hand side for the nonlinear system. 

140 factor : float 

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

142 u0 : dtype_u 

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

144 t : float 

145 Current time. 

146 

147 Returns 

148 ------- 

149 u : dtype_u 

150 Solution. 

151 """ 

152 

153 u = self.dtype_u(u0) 

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

155 b = self.dtype_u(rhs) 

156 

157 self.u_D.t = t 

158 

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

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

161 

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

163 

164 return u 

165 

166 def eval_f(self, u, t): 

167 """ 

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

169 

170 Parameters 

171 ---------- 

172 u : dtype_u 

173 Current values of the numerical solution. 

174 t : float 

175 Current time at which the numerical solution is computed. 

176 

177 Returns 

178 ------- 

179 f : dtype_f 

180 The right-hand side divided into two parts. 

181 """ 

182 

183 f = self.dtype_f(self.V) 

184 

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

186 

187 self.g.t = t 

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

189 return f 

190 

191 def apply_mass_matrix(self, u): 

192 r""" 

193 Routine to apply mass matrix. 

194 

195 Parameters 

196 ---------- 

197 u : dtype_u 

198 Current values of the numerical solution. 

199 

200 Returns 

201 ------- 

202 me : dtype_u 

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

204 """ 

205 

206 me = self.dtype_u(self.V) 

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

208 

209 return me 

210 

211 def u_exact(self, t): 

212 r""" 

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

214 

215 Parameters 

216 ---------- 

217 t : float 

218 Time of the exact solution. 

219 

220 Returns 

221 ------- 

222 me : dtype_u 

223 Exact solution. 

224 """ 

225 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) 

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

227 

228 return me 

229 

230 def fix_residual(self, res): 

231 """ 

232 Applies homogeneous Dirichlet boundary conditions to the residual 

233 

234 Parameters 

235 ---------- 

236 res : dtype_u 

237 Residual 

238 """ 

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

240 return None