Coverage for pySDC/projects/FEM_with_FEniCS/problem_classes/DG_1D_FEniCS.py: 100%

1r""" 

2Discontinuous-Galerkin counterparts of the three CG examples, in mass-matrix form. 

3 

4Each class is its CG parent with the weak form replaced: interior penalty for diffusion, Nitsche for 

5Dirichlet data, local Lax-Friedrichs for the Burgers advection. The Newton loop, mass matrix, initial 

6conditions and solver interface are inherited. 

7 

8The ``penalty`` argument must be pinned to the fine level's value scaled by :math:`h_l / h_0`. The 

9SIPG form depends on the mesh and order it is built on, so rediscretising it per level changes the 

10penalty -- the dominant term -- and the coarse operator stops being :math:`P^T A_F P`. ``setups.py`` 

11builds the pinned ladder. 

12""" 

13 

14import dolfin as df 

15 

16from pySDC.implementations.problem_classes.HeatEquation_1D_FEniCS_matrix_forced import fenics_heat_mass 

17from pySDC.implementations.problem_classes.GrayScott_1D_FEniCS_implicit import fenics_grayscott_mass 

18from pySDC.projects.FEM_with_FEniCS.problem_classes.Burgers_1D_FEniCS import fenics_burgers_mass 

19 

20 

21def sipg(u, v, kappa, mesh, alpha): 

22 r""" 

23 Symmetric interior penalty form for :math:`-\nabla \cdot (\kappa \nabla u)`, natural boundaries. 

24 

25 The exterior facets carry no terms, so this alone imposes homogeneous Neumann conditions. Add 

26 :func:`nitsche` on top for Dirichlet data. 

27 

28 Parameters 

29 ---------- 

30 u : Function, TrialFunction or split component 

31 Argument of the operator. May be nonlinear in ``u``; ``df.derivative`` handles it. 

32 v : TestFunction or split component 

33 Test function. 

34 kappa : Constant or Expression 

35 Diffusion coefficient. 

36 mesh : Mesh 

37 Mesh the form lives on. 

38 alpha : float 

39 Penalty coefficient, entering as :math:`\alpha \kappa / h`. See the module docstring. 

40 

41 Returns 

42 ------- 

43 Form 

44 The bilinear form :math:`a(u, v)`, so a right-hand side uses ``-sipg(...)``. 

45 """ 

46 

47 n, h = df.FacetNormal(mesh), df.CellDiameter(mesh) 

48 alpha = df.Constant(alpha) 

49 

50 return ( 

51 kappa * df.inner(df.grad(u), df.grad(v)) * df.dx 

52 - kappa * df.inner(df.avg(df.grad(u)), df.jump(v, n)) * df.dS 

53 - kappa * df.inner(df.jump(u, n), df.avg(df.grad(v))) * df.dS 

54 + alpha * kappa / df.avg(h) * df.inner(df.jump(u, n), df.jump(v, n)) * df.dS 

55 ) 

56 

57 

58def nitsche(u, v, u_D, kappa, mesh, alpha): 

59 r""" 

60 Weak (Nitsche) Dirichlet terms on the exterior boundary, to be added to :func:`sipg`. 

61 

62 Returns the bilinear part, in :math:`u`, and the linear part carrying the boundary data 

63 :math:`u_D`. They are consistent by construction: for :math:`u = u_D` on the boundary the two 

64 data-dependent terms cancel. 

65 

66 Parameters 

67 ---------- 

68 u, v, kappa, mesh, alpha 

69 As in :func:`sipg`. 

70 u_D : Constant or Expression 

71 Boundary data. 

72 

73 Returns 

74 ------- 

75 a : Form 

76 Bilinear boundary terms, added to the operator. 

77 L : Form 

78 Linear boundary terms, a load vector. 

79 """ 

80 

81 n, h = df.FacetNormal(mesh), df.CellDiameter(mesh) 

82 alpha = df.Constant(alpha) 

83 

84 a = (-kappa * df.dot(df.grad(u), n) * v - kappa * u * df.dot(df.grad(v), n) + alpha * kappa / h * u * v) * df.ds 

85 L = (-kappa * u_D * df.dot(df.grad(v), n) + alpha * kappa / h * u_D * v) * df.ds 

86 

87 return a, L 

88 

89 

90class fenics_heat_dg_mass(fenics_heat_mass): 

91 r""" 

92 Forced 1D heat equation, DG in space, mass-matrix form. 

93 

94 Same problem and same exact solution as :class:`fenics_heat_mass`; the Laplacian becomes SIPG 

95 and the Dirichlet condition :math:`u = c` is imposed by Nitsche terms. 

96 

97 The inherited ``bc`` and ``bc_hom`` survive but have no dofs to act on in a DG space, so 

98 ``solve_system`` and the residual fix are inherited untouched and simply do nothing. 

99 

100 Parameters 

101 ---------- 

102 sigma : float, optional 

103 Interior/boundary penalty constant. Used only when ``penalty`` is not given. 

104 penalty : float, optional 

105 Penalty coefficient, overriding :math:`\sigma p^2`. A multilevel hierarchy has to pin this 

106 across levels, see the module docstring. 

107 family : str, optional 

108 Element family, ``'DG'``. Present so that the level hierarchy can pass it explicitly. 

109 **kwargs 

110 Forwarded to :class:`fenics_heat_mass`. 

111 

112 Attributes 

113 ---------- 

114 b_bc : GenericVector 

115 Assembled Nitsche load vector carrying the boundary data. 

116 """ 

117 

118 def __init__(self, sigma=10.0, penalty=None, family='DG', **kwargs): 

119 """Initialization routine""" 

120 

121 super().__init__(family=family, **kwargs) 

122 self._makeAttributeAndRegister('sigma', 'penalty', localVars=locals(), readOnly=True) 

123 alpha = sigma * max(self.order, 1) ** 2 if penalty is None else penalty 

124 

125 # no strong boundary conditions to fix the residual against 

126 self.fix_bc_for_residual = False 

127 

128 mesh = self.V.mesh() 

129 u, v = df.TrialFunction(self.V), df.TestFunction(self.V) 

130 kappa = df.Constant(self.nu) 

131 a_bc, L_bc = nitsche(u, v, df.Constant(self.c), kappa, mesh, alpha) 

132 

133 # K is the negated operator, matching the parent: M u' = K u + b_bc + M g 

134 self.K = df.assemble(-(sipg(u, v, kappa, mesh, alpha) + a_bc)) 

135 self.b_bc = df.assemble(L_bc) 

136 

137 def eval_f(self, u, t): 

138 """ 

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

140 

141 The Nitsche load rides in the explicit part rather than in the implicit one, which keeps 

142 ``solve_system`` identical to the parent's. It is constant in time here, and every 

143 quadrature involved integrates constants exactly, so the split costs nothing: a constant 

144 does not enter the iteration matrix. 

145 

146 Parameters 

147 ---------- 

148 u : dtype_u 

149 Current values of the numerical solution. 

150 t : float 

151 Current time at which the numerical solution is computed. 

152 

153 Returns 

154 ------- 

155 f : dtype_f 

156 The right-hand side divided into two parts. 

157 """ 

158 

159 f = super().eval_f(u, t) 

160 f.expl.values.vector().axpy(1.0, self.b_bc) 

161 

162 return f 

163 

164 

165class fenics_burgers_dg_mass(fenics_burgers_mass): 

166 r""" 

167 Viscous Burgers in 1D, DG in space, mass-matrix form. 

168 

169 The advection is taken in conservative form, :math:`(u^2/2)_x`, with a local Lax-Friedrichs 

170 numerical flux on the interior facets. The exterior facets use the interior trace, which leaves 

171 the boundary transmissive -- the same "no condition imposed" situation as the non-conservative 

172 CG form, and harmless for the cosine initial condition over the times run here. 

173 

174 The viscous term is SIPG with natural boundaries. 

175 

176 Parameters 

177 ---------- 

178 sigma : float, optional 

179 Interior penalty constant. Used only when ``penalty`` is not given. 

180 penalty : float, optional 

181 Penalty coefficient, overriding :math:`\sigma p^2`. A multilevel hierarchy has to pin this 

182 across levels, see the module docstring. 

183 family : str, optional 

184 Element family, ``'DG'``. 

185 **kwargs 

186 Forwarded to :class:`fenics_burgers_mass`. 

187 """ 

188 

189 def __init__(self, sigma=10.0, penalty=None, family='DG', **kwargs): 

190 """Initialization routine""" 

191 

192 super().__init__(family=family, **kwargs) 

193 self._makeAttributeAndRegister('sigma', 'penalty', localVars=locals(), readOnly=True) 

194 alpha = sigma * max(self.order, 1) ** 2 if penalty is None else penalty 

195 

196 mesh = self.V.mesh() 

197 n = df.FacetNormal(mesh) 

198 v = df.TestFunction(self.V) 

199 w = self.w 

200 

201 flux = df.as_vector([w**2 / 2]) 

202 left, right = abs(w('+')), abs(w('-')) 

203 speed = df.conditional(df.gt(left, right), left, right) 

204 flux_hat = df.avg(flux) + 0.5 * speed * df.jump(w, n) 

205 

206 self.F = ( 

207 df.dot(flux, df.grad(v)) * df.dx 

208 - df.dot(flux_hat, df.jump(v, n)) * df.dS 

209 - df.dot(flux, n) * v * df.ds 

210 - sipg(w, v, df.Constant(self.nu), mesh, alpha) 

211 ) 

212 

213 

214class fenics_grayscott_dg_mass(fenics_grayscott_mass): 

215 r""" 

216 Gray-Scott in 1D, DG in space, mass-matrix form. 

217 

218 Both diffusion terms become SIPG with natural boundaries, matching the CG version's Neumann 

219 conditions. The reaction terms are cell-local and carry over verbatim. 

220 

221 Parameters 

222 ---------- 

223 sigma : float, optional 

224 Interior penalty constant. Used only when ``penalty`` is not given. 

225 penalty : float, optional 

226 Penalty coefficient, overriding :math:`\sigma p^2`. A multilevel hierarchy has to pin this 

227 across levels, see the module docstring. 

228 family : str, optional 

229 Element family, ``'DG'``. 

230 **kwargs 

231 Forwarded to :class:`fenics_grayscott_mass`. 

232 """ 

233 

234 def __init__(self, sigma=10.0, penalty=None, family='DG', **kwargs): 

235 """Initialization routine""" 

236 

237 super().__init__(family=family, **kwargs) 

238 self._makeAttributeAndRegister('sigma', 'penalty', localVars=locals(), readOnly=True) 

239 alpha = sigma * max(self.order, 1) ** 2 if penalty is None else penalty 

240 

241 mesh = self.V.mesh() 

242 q1, q2 = df.TestFunctions(self.V) 

243 w1, w2 = self.w1, self.w2 

244 

245 self.F1 = (-w1 * (w2**2) * q1 + self.A * (1 - w1) * q1) * df.dx - sipg( 

246 w1, q1, df.Constant(self.Du), mesh, alpha 

247 ) 

248 self.F2 = (w1 * (w2**2) * q2 - self.B * w2 * q2) * df.dx - sipg(w2, q2, df.Constant(self.Dv), mesh, alpha) 

249 self.F = self.F1 + self.F2