Coverage for pySDC/implementations/transfer_classes/TransferFenicsMesh.py: 100%

1import dolfin as df 

2import numpy as np 

3import scipy.sparse as sp 

4import scipy.sparse.linalg as spla 

5 

6from pySDC.core.errors import TransferError 

7from pySDC.core.space_transfer import SpaceTransfer 

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

9 

10 

11class mesh_to_mesh_fenics(SpaceTransfer): 

12 """ 

13 This implementation can restrict and prolong between fenics meshes 

14 """ 

15 

16 def __init__(self, fine_prob, coarse_prob, params): 

17 super().__init__(fine_prob, coarse_prob, params) 

18 self._Pmat = None 

19 self._l2 = None 

20 

21 @property 

22 def Pmat(self): 

23 """ 

24 Prolongation matrix P (coarse -> fine), assembled once on first use. 

25 

26 P is the inclusion of the coarse space in the fine one: column j holds the coarse basis 

27 function phi_j expanded in the fine basis. It is built cell by cell -- for every fine cell 

28 the enclosing coarse cell is found from the fine cell's midpoint, which is never on a coarse 

29 facet, and the coarse basis is evaluated in *that* cell. 

30 

31 Not df.interpolate: across meshes that is point evaluation, and a fine dof on a coarse facet 

32 has two coarse values there, so a discontinuous coarse function comes back continuised. For 

33 continuous spaces the two agree to machine precision. 

34 

35 Costs one basis evaluation per fine dof; dolfin 2019.1.0's 

36 PETScDMCollection.create_transfer_matrix segfaults. 

37 """ 

38 if self._Pmat is None: 

39 Vc, Vf = self.coarse_prob.init, self.fine_prob.init 

40 tree = Vc.mesh().bounding_box_tree() 

41 element, dofmap_c = Vc.element(), Vc.dofmap() 

42 element_f = Vf.element() 

43 

44 # which component of a mixed space each fine dof belongs to; scalar spaces are all zero 

45 ncomp = max(Vf.num_sub_spaces(), 1) 

46 component = np.zeros(Vf.dim(), dtype=int) 

47 for k in range(Vf.num_sub_spaces()): 

48 component[Vf.sub(k).dofmap().dofs()] = k 

49 

50 rows, cols, vals = [], [], [] 

51 seen = set() 

52 for cell_f in df.cells(Vf.mesh()): 

53 cell_c = df.Cell(Vc.mesh(), tree.compute_first_entity_collision(cell_f.midpoint())) 

54 coords, orientation = cell_c.get_vertex_coordinates(), cell_c.orientation() 

55 dofs_c = dofmap_c.cell_dofs(cell_c.index()) 

56 # per cell, not from tabulate_dof_coordinates: that reports one point per master 

57 # dof, which on a periodic space is the wrong side of the domain for half the cells 

58 # and evaluates the coarse basis outside the cell 

59 x_local = element_f.tabulate_dof_coordinates(cell_f) 

60 for k, dof_f in enumerate(Vf.dofmap().cell_dofs(cell_f.index())): 

61 # a continuous space shares dofs between cells; the second visit is redundant 

62 if dof_f in seen: 

63 continue 

64 seen.add(dof_f) 

65 basis = np.asarray(element.evaluate_basis_all(x_local[k], coords, orientation)) 

66 col = basis.reshape(-1, ncomp)[:, component[dof_f]] 

67 nz = np.nonzero(np.abs(col) > 1e-13)[0] 

68 rows.extend([dof_f] * len(nz)) 

69 cols.extend(dofs_c[nz]) 

70 vals.extend(col[nz]) 

71 self._Pmat = sp.csr_matrix((vals, (rows, cols)), shape=(Vf.dim(), Vc.dim())) 

72 return self._Pmat 

73 

74 def restrict_dual(self, F): 

75 """ 

76 Variational restriction P^T, for quantities living in the dual space. 

77 

78 The FAS tau of the mass formulation is a load vector, not a nodal function, so it has to 

79 be restricted with P^T rather than by interpolation. For nested Lagrange spaces 

80 phi_i^coarse = sum_j P_ji phi_j^fine, which makes P^T exact: it reproduces the coarse 

81 load vector of the same functional. Interpolating instead reads the load vector as if it 

82 were a function and is wrong by roughly 2^dim. 

83 

84 Args: 

85 F: the fine level data 

86 """ 

87 if not isinstance(F, fenics_mesh): 

88 # tau and u0 are the only dual quantities, and both are meshes 

89 raise TransferError('Unknown type of fine data, got %s' % type(F)) 

90 

91 u_coarse = fenics_mesh(self.coarse_prob.init) 

92 u_coarse.values.vector()[:] = self.Pmat.T.dot(F.values.vector()[:]) 

93 return u_coarse 

94 

95 @property 

96 def l2_pieces(self): 

97 """ 

98 Prefactorised coarse mass matrix, :math:`P^T`, and the fine mass matrix, assembled once. 

99 

100 Factorised up front: it is solved against once per node per sweep. 

101 """ 

102 if self._l2 is None: 

103 

104 def mass(V): 

105 # straight out of PETSc's CSR; `.array()` would densify, and these are ~0.5% dense 

106 u, v = df.TrialFunction(V), df.TestFunction(V) 

107 mat = df.as_backend_type(df.assemble(df.inner(u, v) * df.dx)).mat() 

108 indptr, indices, data = mat.getValuesCSR() 

109 return sp.csr_matrix((data, indices, indptr), shape=mat.getSize()) 

110 

111 self._l2 = (spla.factorized(mass(self.coarse_prob.init).tocsc()), self.Pmat.T, mass(self.fine_prob.init)) 

112 return self._l2 

113 

114 def _project_one(self, values): 

115 """L2-project one fine ``dolfin.Function`` onto the coarse space.""" 

116 solve_coarse, P_transpose, mass_fine = self.l2_pieces 

117 coarse = df.Function(self.coarse_prob.init) 

118 coarse.vector()[:] = solve_coarse(P_transpose @ (mass_fine @ values.vector()[:])) 

119 return coarse 

120 

121 def project(self, F): 

122 """ 

123 Restriction of a SOLUTION, by L2 projection: :math:`M_c^{-1} P^T M_f`. 

124 

125 Not point sampling, which is cheaper to reach for but is not well defined on a DG space and 

126 costs PFASST iterations even on a CG one -- :math:`R_u` cancels out of the linear FAS 

127 iteration but not across step boundaries. See ``projects/FEM_with_FEniCS`` for the numbers. 

128 

129 Dual quantities go through :meth:`restrict_dual` instead. 

130 

131 Args: 

132 F: the fine level data 

133 """ 

134 if isinstance(F, fenics_mesh): 

135 return fenics_mesh(self._project_one(F.values)) 

136 elif isinstance(F, rhs_fenics_mesh): 

137 u_coarse = rhs_fenics_mesh(self.coarse_prob.init) 

138 u_coarse.impl.values = self._project_one(F.impl.values) 

139 u_coarse.expl.values = self._project_one(F.expl.values) 

140 return u_coarse 

141 raise TransferError('Unknown type of fine data, got %s' % type(F)) 

142 

143 def restrict(self, F): 

144 """ 

145 Restriction implementation 

146 

147 Args: 

148 F: the fine level data 

149 """ 

150 if isinstance(F, fenics_mesh): 

151 u_coarse = fenics_mesh(df.interpolate(F.values, self.coarse_prob.init)) 

152 elif isinstance(F, rhs_fenics_mesh): 

153 u_coarse = rhs_fenics_mesh(self.coarse_prob.init) 

154 u_coarse.impl.values = df.interpolate(F.impl.values, self.coarse_prob.init) 

155 u_coarse.expl.values = df.interpolate(F.expl.values, self.coarse_prob.init) 

156 else: 

157 raise TransferError('Unknown type of fine data, got %s' % type(F)) 

158 

159 return u_coarse 

160 

161 def prolong(self, G): 

162 """ 

163 Prolongation implementation, the exact inclusion via P. 

164 

165 Not df.interpolate: see the Pmat docstring for why that silently continuises a 

166 discontinuous coarse function. 

167 

168 Args: 

169 G: the coarse level data 

170 """ 

171 P = self.Pmat 

172 if isinstance(G, fenics_mesh): 

173 u_fine = fenics_mesh(self.fine_prob.init) 

174 u_fine.values.vector()[:] = P.dot(G.values.vector()[:]) 

175 elif isinstance(G, rhs_fenics_mesh): 

176 u_fine = rhs_fenics_mesh(self.fine_prob.init) 

177 u_fine.impl.values.vector()[:] = P.dot(G.impl.values.vector()[:]) 

178 u_fine.expl.values.vector()[:] = P.dot(G.expl.values.vector()[:]) 

179 else: 

180 raise TransferError('Unknown type of coarse data, got %s' % type(G)) 

181 

182 return u_fine