Coverage for pySDC/projects/FEM_with_FEniCS/problem_classes/Burgers_1D_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 

8 

9 

10class fenics_burgers_mass(Problem): 

11 r""" 

12 Viscous Burgers equation in one dimension, in mass-matrix form. 

13 

14 .. math:: 

15 \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} 

16 = \nu \frac{\partial^2 u}{\partial x^2} 

17 

18 on :math:`[0, 1]` with natural (homogeneous Neumann) boundary conditions, discretised with 

19 ``FEniCS``. In weak form, with :math:`M` the mass matrix, 

20 

21 .. math:: 

22 M \vec{u}' = F(\vec{u}), \qquad 

23 F(\vec{u}) = -\left(u u_x, v\right) - \nu \left(u_x, v_x\right). 

24 

25 ``eval_f`` returns the assembled :math:`F`, a load vector, and ``solve_system`` takes a 

26 right-hand side that is already in the dual space. The mass matrix is never inverted. 

27 

28 The nonlinearity is quadratic *advection* rather than a reaction term, which is what makes this 

29 example worth having next to Gray-Scott: advection is where the order of the prolongation stops 

30 being a tuning knob and becomes a stability requirement. 

31 

32 Parameters 

33 ---------- 

34 c_nvars : int, optional 

35 Number of cells on the base mesh. 

36 family : str, optional 

37 Finite element family. 

38 order : int, optional 

39 Element order. High order is the point here: the solution is smooth for moderate ``nu``, and 

40 mass lumping is not available at this order anyway. 

41 refinements : int, optional 

42 Number of mesh refinements, one entry per level for a multilevel run. 

43 nu : float, optional 

44 Viscosity. Keep it moderate: as it shrinks the solution develops a steep front, high-order 

45 CG starts to oscillate without stabilisation, and the coarse level stops being able to carry 

46 the solution. 

47 newton_tol, newton_rtol : float, optional 

48 Absolute and relative tolerance of the node-local Newton solve. Both have to be tighter than 

49 the SDC residual tolerance or the node-local solve sets the accuracy floor. 

50 newton_maxiter : int, optional 

51 Maximum number of node-local Newton iterations. 

52 

53 Attributes 

54 ---------- 

55 V : FunctionSpace 

56 Function space of the problem. 

57 M : scalar 

58 Mass matrix. 

59 w : Function 

60 Working function carrying the current Newton iterate; ``F`` is written in terms of it. 

61 F : Form 

62 Weak form of the right-hand side. 

63 """ 

64 

65 dtype_u = fenics_mesh 

66 dtype_f = fenics_mesh 

67 

68 def __init__( 

69 self, 

70 c_nvars=128, 

71 family='CG', 

72 order=4, 

73 refinements=0, 

74 nu=0.02, 

75 newton_tol=1e-12, 

76 newton_rtol=1e-11, 

77 newton_maxiter=100, 

78 ): 

79 """Initialization routine""" 

80 

81 warning_level = getattr(df, 'WARNING', None) 

82 if warning_level is None and hasattr(df, 'LogLevel'): 

83 warning_level = df.LogLevel.WARNING 

84 if warning_level is not None: 

85 df.set_log_level(warning_level) 

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

87 

88 df.parameters['form_compiler']['optimize'] = True 

89 df.parameters['form_compiler']['cpp_optimize'] = True 

90 

91 mesh = df.UnitIntervalMesh(c_nvars) 

92 num_refinements = refinements if isinstance(refinements, int) else sum(refinements) 

93 for _ in range(num_refinements): 

94 mesh = df.refine(mesh) 

95 

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

97 

98 super().__init__(self.V) 

99 self._makeAttributeAndRegister( 

100 'c_nvars', 'family', 'order', 'refinements', 'nu', localVars=locals(), readOnly=True 

101 ) 

102 self._makeAttributeAndRegister( 

103 'newton_tol', 'newton_rtol', 'newton_maxiter', localVars=locals(), readOnly=False 

104 ) 

105 

106 self.w = df.Function(self.V) 

107 v = df.TestFunction(self.V) 

108 self.F = (-self.w * self.w.dx(0) * v - self.nu * self.w.dx(0) * v.dx(0)) * df.dx 

109 

110 u, q = df.TrialFunction(self.V), df.TestFunction(self.V) 

111 self.M = df.assemble(u * q * df.dx) 

112 

113 def apply_mass_matrix(self, u): 

114 r""" 

115 Apply the mass matrix, :math:`M \vec{u}`. 

116 

117 Parameters 

118 ---------- 

119 u : dtype_u 

120 Current values of the numerical solution. 

121 

122 Returns 

123 ------- 

124 me : dtype_u 

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

126 """ 

127 

128 me = self.dtype_u(self.V) 

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

130 

131 return me 

132 

133 def eval_f(self, u, t): 

134 r""" 

135 Evaluate the right-hand side in weak (dual) form, without inverting the mass matrix. 

136 

137 Parameters 

138 ---------- 

139 u : dtype_u 

140 Current values of the numerical solution. 

141 t : float 

142 Current time. 

143 

144 Returns 

145 ------- 

146 f : dtype_f 

147 The assembled load vector :math:`F(\vec{u})`. 

148 """ 

149 

150 f = self.dtype_f(self.V) 

151 self.w.assign(u.values) 

152 f.values = df.Function(self.V, df.assemble(self.F)) 

153 

154 return f 

155 

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

157 r""" 

158 Solve :math:`M \vec{u} - factor \cdot F(\vec{u}) = \vec{rhs}` by Newton, rhs already dual. 

159 

160 Parameters 

161 ---------- 

162 rhs : dtype_f 

163 Right-hand side, in the dual space. 

164 factor : float 

165 Node-to-node stepsize. 

166 u0 : dtype_u 

167 Initial guess for the Newton iteration. 

168 t : float 

169 Current time. 

170 

171 Returns 

172 ------- 

173 sol : dtype_u 

174 Solution. 

175 """ 

176 

177 sol = self.dtype_u(u0) 

178 self.w.assign(sol.values) 

179 

180 du = df.TrialFunction(self.V) 

181 Jform = df.derivative(self.F, self.w, du) 

182 

183 res, delta, res0 = df.Function(self.V), df.Function(self.V), None 

184 

185 for _ in range(self.newton_maxiter): 

186 self.M.mult(self.w.vector(), res.vector()) 

187 res.vector().axpy(-factor, df.assemble(self.F)) 

188 res.vector().axpy(-1.0, rhs.values.vector()) 

189 

190 norm = res.vector().norm('l2') 

191 res0 = norm if res0 is None else res0 

192 if norm < self.newton_tol or norm < self.newton_rtol * res0: 

193 break 

194 

195 df.solve(self.M - factor * df.assemble(Jform), delta.vector(), res.vector()) 

196 self.w.vector().axpy(-1.0, delta.vector()) 

197 

198 sol.values.assign(self.w) 

199 

200 return sol 

201 

202 def u_exact(self, t): 

203 r""" 

204 Initial condition :math:`u(x, 0) = \cos(2 \pi x)`, only available at :math:`t = 0`. 

205 

206 It has vanishing derivative at both ends, so it is consistent with the natural boundary 

207 conditions, and it is smooth, so the coarse levels can carry it. 

208 

209 Parameters 

210 ---------- 

211 t : float 

212 Time of the exact solution. 

213 

214 Returns 

215 ------- 

216 me : dtype_u 

217 Initial condition. 

218 """ 

219 

220 assert t == 0, 'ERROR: u_exact only valid for t=0' 

221 

222 expr = df.Expression('cos(2*pi*x[0])', degree=max(1, self.order), pi=np.pi) 

223 

224 return self.dtype_u(df.interpolate(expr, self.V))