Coverage for pySDC/projects/AllenCahn_Bayreuth/AllenCahn_monitor.py: 87%

1import numpy as np 

2from mpi4py import MPI 

3from pySDC.core.Hooks import hooks 

4 

5 

6class monitor(hooks): 

7 def __init__(self): 

8 """ 

9 Initialization of Allen-Cahn monitoring 

10 """ 

11 super(monitor, self).__init__() 

12 

13 self.init_radius = None 

14 self.init_vol = None 

15 self.ndim = None 

16 self.corr_rad = None 

17 self.corr_vol = None 

18 

19 self.comm = None 

20 

21 def pre_run(self, step, level_number): 

22 """ 

23 Overwrite standard pre run hook 

24 

25 Args: 

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

27 level_number (int): the current level number 

28 """ 

29 super(monitor, self).pre_run(step, level_number) 

30 L = step.levels[0] 

31 

32 # get space-communicator and data 

33 self.comm = L.prob.comm 

34 

35 # get real space values 

36 if L.prob.spectral: 

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

38 else: 

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

40 self.ndim = len(tmp.shape) 

41 

42 # compute numerical radius and volume 

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

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

45 if self.comm is not None: 

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

47 else: 

48 c_global = c_local 

49 if self.ndim == 3: 

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

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

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

53 elif self.ndim == 2: 

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

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

56 self.init_vol = np.pi * L.prob.radius**2 

57 else: 

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

59 

60 self.init_radius = L.prob.radius 

61 self.corr_rad = self.init_radius / radius 

62 self.corr_vol = self.init_vol / vol 

63 radius *= self.corr_rad 

64 vol *= self.corr_vol 

65 

66 # write to stats 

67 if L.time == 0.0: 

68 self.add_to_stats( 

69 process=step.status.slot, 

70 time=L.time, 

71 level=-1, 

72 iter=step.status.iter, 

73 sweep=L.status.sweep, 

74 type='computed_radius', 

75 value=radius, 

76 ) 

77 self.add_to_stats( 

78 process=step.status.slot, 

79 time=L.time, 

80 level=-1, 

81 iter=step.status.iter, 

82 sweep=L.status.sweep, 

83 type='exact_radius', 

84 value=self.init_radius, 

85 ) 

86 self.add_to_stats( 

87 process=step.status.slot, 

88 time=L.time, 

89 level=-1, 

90 iter=step.status.iter, 

91 sweep=L.status.sweep, 

92 type='computed_volume', 

93 value=vol, 

94 ) 

95 self.add_to_stats( 

96 process=step.status.slot, 

97 time=L.time, 

98 level=-1, 

99 iter=step.status.iter, 

100 sweep=L.status.sweep, 

101 type='exact_volume', 

102 value=self.init_vol, 

103 ) 

104 

105 def post_step(self, step, level_number): 

106 """ 

107 Overwrite standard post step hook 

108 

109 Args: 

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

111 level_number (int): the current level number 

112 """ 

113 super(monitor, self).post_step(step, level_number) 

114 

115 # some abbreviations 

116 L = step.levels[0] 

117 

118 # get real space values 

119 if L.prob.spectral: 

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

121 else: 

122 tmp = L.uend[:] 

123 

124 # compute numerical radius and volume 

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

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

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

128 if self.comm is not None: 

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

130 else: 

131 c_global = c_local 

132 

133 if self.ndim == 3: 

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

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

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

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

138 elif self.ndim == 2: 

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

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

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

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

143 else: 

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

145 

146 radius *= self.corr_rad 

147 vol *= self.corr_vol 

148 

149 # write to stats 

150 self.add_to_stats( 

151 process=step.status.slot, 

152 time=L.time + L.dt, 

153 level=-1, 

154 iter=step.status.iter, 

155 sweep=L.status.sweep, 

156 type='computed_radius', 

157 value=radius, 

158 ) 

159 self.add_to_stats( 

160 process=step.status.slot, 

161 time=L.time + L.dt, 

162 level=-1, 

163 iter=step.status.iter, 

164 sweep=L.status.sweep, 

165 type='exact_radius', 

166 value=exact_radius, 

167 ) 

168 self.add_to_stats( 

169 process=step.status.slot, 

170 time=L.time + L.dt, 

171 level=-1, 

172 iter=step.status.iter, 

173 sweep=L.status.sweep, 

174 type='computed_volume', 

175 value=vol, 

176 ) 

177 self.add_to_stats( 

178 process=step.status.slot, 

179 time=L.time + L.dt, 

180 level=-1, 

181 iter=step.status.iter, 

182 sweep=L.status.sweep, 

183 type='exact_volume', 

184 value=exact_vol, 

185 )