Coverage for pySDC/projects/AllenCahn_Bayreuth/AllenCahn_monitor_and_dump.py: 0%

1import numpy as np 

2import json 

3from mpi4py import MPI 

4from pySDC.core.Hooks import hooks 

5 

6 

7class monitor_and_dump(hooks): 

8 def __init__(self): 

9 """ 

10 Initialization of Allen-Cahn monitoring 

11 """ 

12 super(monitor_and_dump, self).__init__() 

13 

14 self.init_radius = None 

15 self.init_vol = None 

16 self.ndim = None 

17 self.corr_rad = None 

18 self.corr_vol = None 

19 

20 self.comm = None 

21 self.rank = None 

22 self.size = None 

23 self.amode = MPI.MODE_WRONLY | MPI.MODE_CREATE 

24 self.time_step = None 

25 

26 def pre_run(self, step, level_number): 

27 """ 

28 Overwrite standard pre run hook 

29 

30 Args: 

31 step (pySDC.Step.step): the current step 

32 level_number (int): the current level number 

33 """ 

34 super(monitor_and_dump, self).pre_run(step, level_number) 

35 L = step.levels[0] 

36 

37 # get space-communicator and data 

38 self.comm = L.prob.params.comm 

39 if self.comm is not None: 

40 self.rank = self.comm.Get_rank() 

41 self.size = self.comm.Get_size() 

42 else: 

43 self.rank = 0 

44 self.size = 1 

45 

46 # get real space values 

47 if L.prob.params.spectral: 

48 tmp = L.prob.fft.backward(L.u[0]) 

49 else: 

50 tmp = L.u[0][:] 

51 

52 self.ndim = len(tmp.shape) 

53 

54 # compute numerical radius and volume 

55 # c_local = np.count_nonzero(tmp >= 0.5) 

56 c_local = float(tmp[:].sum()) 

57 if self.comm is not None: 

58 c_global = self.comm.allreduce(sendobj=c_local, op=MPI.SUM) 

59 else: 

60 c_global = c_local 

61 if self.ndim == 3: 

62 vol = c_global * L.prob.dx**3 

63 radius = (vol / (np.pi * 4.0 / 3.0)) ** (1.0 / 3.0) 

64 self.init_vol = np.pi * 4.0 / 3.0 * L.prob.params.radius**3 

65 elif self.ndim == 2: 

66 vol = c_global * L.prob.dx**2 

67 radius = np.sqrt(vol / np.pi) 

68 self.init_vol = np.pi * L.prob.params.radius**2 

69 else: 

70 raise NotImplementedError('Can use this only for 2 or 3D problems') 

71 

72 self.init_radius = L.prob.params.radius 

73 self.corr_rad = self.init_radius / radius 

74 self.corr_vol = self.init_vol / vol 

75 radius *= self.corr_rad 

76 vol *= self.corr_vol 

77 

78 # write to stats 

79 if L.time == 0.0: 

80 self.add_to_stats( 

81 process=step.status.slot, 

82 time=L.time, 

83 level=-1, 

84 iter=step.status.iter, 

85 sweep=L.status.sweep, 

86 type='computed_radius', 

87 value=radius, 

88 ) 

89 self.add_to_stats( 

90 process=step.status.slot, 

91 time=L.time, 

92 level=-1, 

93 iter=step.status.iter, 

94 sweep=L.status.sweep, 

95 type='exact_radius', 

96 value=self.init_radius, 

97 ) 

98 self.add_to_stats( 

99 process=step.status.slot, 

100 time=L.time, 

101 level=-1, 

102 iter=step.status.iter, 

103 sweep=L.status.sweep, 

104 type='computed_volume', 

105 value=vol, 

106 ) 

107 self.add_to_stats( 

108 process=step.status.slot, 

109 time=L.time, 

110 level=-1, 

111 iter=step.status.iter, 

112 sweep=L.status.sweep, 

113 type='exact_volume', 

114 value=self.init_vol, 

115 ) 

116 

117 # compute local offset for I/O 

118 nbytes_local = tmp.nbytes 

119 if self.comm is not None: 

120 nbytes_global = self.comm.allgather(nbytes_local) 

121 else: 

122 nbytes_global = [nbytes_local] 

123 local_offset = sum(nbytes_global[: self.rank]) 

124 

125 # dump initial data 

126 fname = f"./data/{L.prob.params.name}_{0:08d}" 

127 fh = MPI.File.Open(self.comm, fname + ".dat", self.amode) 

128 fh.Write_at_all(local_offset, tmp) 

129 fh.Close() 

130 

131 # write json description 

132 if self.rank == 0 and step.status.slot == 0: 

133 json_obj = dict() 

134 json_obj['type'] = 'dataset' 

135 json_obj['datatype'] = str(tmp.dtype) 

136 json_obj['endian'] = str(tmp.dtype.byteorder) 

137 json_obj['time'] = L.time 

138 json_obj['space_comm_size'] = self.size 

139 json_obj['time_comm_size'] = step.status.time_size 

140 json_obj['shape'] = L.prob.params.nvars 

141 json_obj['elementsize'] = tmp.dtype.itemsize 

142 

143 with open(fname + '.json', 'w') as fp: 

144 json.dump(json_obj, fp) 

145 

146 # set step count 

147 self.time_step = 1 

148 

149 def post_step(self, step, level_number): 

150 """ 

151 Overwrite standard post step hook 

152 

153 Args: 

154 step (pySDC.Step.step): the current step 

155 level_number (int): the current level number 

156 """ 

157 super(monitor_and_dump, self).post_step(step, level_number) 

158 

159 # some abbreviations 

160 L = step.levels[0] 

161 

162 # get real space values 

163 if L.prob.params.spectral: 

164 tmp = L.prob.fft.backward(L.uend) 

165 else: 

166 tmp = L.uend[:] 

167 

168 # compute numerical radius and volume 

169 # c_local = np.count_nonzero(tmp >= 0.5) 

170 # c_local = float(tmp[tmp > 2 * L.prob.params.eps].sum()) 

171 c_local = float(tmp[:].sum()) 

172 if self.comm is not None: 

173 c_global = self.comm.allreduce(sendobj=c_local, op=MPI.SUM) 

174 else: 

175 c_global = c_local 

176 

177 if self.ndim == 3: 

178 vol = c_global * L.prob.dx**3 

179 radius = (vol / (np.pi * 4.0 / 3.0)) ** (1.0 / 3.0) 

180 exact_vol = np.pi * 4.0 / 3.0 * (max(self.init_radius**2 - 4.0 * (L.time + L.dt), 0)) ** (3.0 / 2.0) 

181 exact_radius = (exact_vol / (np.pi * 4.0 / 3.0)) ** (1.0 / 3.0) 

182 elif self.ndim == 2: 

183 vol = c_global * L.prob.dx**2 

184 radius = np.sqrt(vol / np.pi) 

185 exact_vol = np.pi * max(self.init_radius**2 - 2.0 * (L.time + L.dt), 0) 

186 exact_radius = np.sqrt(exact_vol / np.pi) 

187 else: 

188 raise NotImplementedError('Can use this only for 2 or 3D problems') 

189 

190 radius *= self.corr_rad 

191 vol *= self.corr_vol 

192 

193 # write to stats 

194 self.add_to_stats( 

195 process=step.status.slot, 

196 time=L.time + L.dt, 

197 level=-1, 

198 iter=step.status.iter, 

199 sweep=L.status.sweep, 

200 type='computed_radius', 

201 value=radius, 

202 ) 

203 self.add_to_stats( 

204 process=step.status.slot, 

205 time=L.time + L.dt, 

206 level=-1, 

207 iter=step.status.iter, 

208 sweep=L.status.sweep, 

209 type='exact_radius', 

210 value=exact_radius, 

211 ) 

212 self.add_to_stats( 

213 process=step.status.slot, 

214 time=L.time + L.dt, 

215 level=-1, 

216 iter=step.status.iter, 

217 sweep=L.status.sweep, 

218 type='computed_volume', 

219 value=vol, 

220 ) 

221 self.add_to_stats( 

222 process=step.status.slot, 

223 time=L.time + L.dt, 

224 level=-1, 

225 iter=step.status.iter, 

226 sweep=L.status.sweep, 

227 type='exact_volume', 

228 value=exact_vol, 

229 ) 

230 

231 # compute local offset for I/O 

232 nbytes_local = tmp.nbytes 

233 if self.comm is not None: 

234 nbytes_global = self.comm.allgather(nbytes_local) 

235 else: 

236 nbytes_global = [nbytes_local] 

237 local_offset = sum(nbytes_global[: self.rank]) 

238 

239 # dump initial data 

240 fname = f"./data/{L.prob.params.name}_{self.time_step + step.status.slot:08d}" 

241 fh = MPI.File.Open(self.comm, fname + ".dat", self.amode) 

242 fh.Write_at_all(local_offset, tmp) 

243 fh.Close() 

244 

245 # write json description 

246 if self.rank == 0: 

247 json_obj = dict() 

248 json_obj['type'] = 'dataset' 

249 json_obj['datatype'] = str(tmp.dtype) 

250 json_obj['endian'] = str(tmp.dtype.byteorder) 

251 json_obj['time'] = L.time + L.dt 

252 json_obj['space_comm_size'] = self.size 

253 json_obj['time_comm_size'] = step.status.time_size 

254 json_obj['shape'] = L.prob.params.nvars 

255 json_obj['elementsize'] = tmp.dtype.itemsize 

256 

257 with open(fname + '.json', 'w') as fp: 

258 json.dump(json_obj, fp) 

259 

260 # update step count 

261 self.time_step += step.status.time_size