# Coverage for pySDC/implementations/problem_classes/AllenCahn_MPIFFT.py: 73%

## 92 statements

, created at 2024-09-19 09:13 +0000

1import numpy as np

2from mpi4py import MPI

4from pySDC.implementations.problem_classes.generic_MPIFFT_Laplacian import IMEX_Laplacian_MPIFFT

5from mpi4py_fft import newDistArray

8class allencahn_imex(IMEX_Laplacian_MPIFFT):

9 r"""

10 Example implementing the :math:2-dimensional Allen-Cahn equation with periodic boundary conditions :math:u \in [0, 1]^2

12 .. math::

13 \frac{\partial u}{\partial t} = \Delta u - \frac{2}{\varepsilon^2} u (1 - u) (1 - 2u)

14 - 6 d_w u (1 - u)

16 on a spatial domain :math:[-\frac{L}{2}, \frac{L}{2}]^2, driving force :math:d_w, and :math:N=2,3. Different initial

17 conditions can be used, for example, circles of the form

19 .. math::

20 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{(x_i-0.5)^2 + (y_j-0.5)^2}}{\sqrt{2}\varepsilon}\right),

22 for :math:i, j=0,..,N-1, where :math:N is the number of spatial grid points. For time-stepping, the problem is treated

23 *semi-implicitly*, i.e., the linear part is solved with Fast-Fourier Transform (FFT) and the nonlinear part in the right-hand

24 side will be treated explicitly using mpi4py-fft [1]_ to solve them.

26 Parameters

27 ----------

28 nvars : List of int tuples, optional

29 Number of unknowns in the problem, e.g. nvars=(128, 128).

30 eps : float, optional

31 Scaling parameter :math:\varepsilon.

34 spectral : bool, optional

35 Indicates if spectral initial condition is used.

36 dw : float, optional

37 Driving force.

38 L : float, optional

39 Denotes the period of the function to be approximated for the Fourier transform.

40 init_type : str, optional

41 Initialises type of initial state.

42 comm : bool, optional

43 Communicator for parallelization.

45 Attributes

46 ----------

47 fft : fft object

48 Object for FFT.

49 X : np.ogrid

50 Grid coordinates in real space.

51 K2 : np.1darray

52 Laplace operator in spectral space.

53 dx : float

54 Mesh width in x direction.

55 dy : float

56 Mesh width in y direction.

58 References

59 ----------

61 """

63 def __init__(

64 self,

65 eps=0.04,

67 dw=0.0,

68 init_type='circle',

69 **kwargs,

70 ):

71 kwargs['L'] = kwargs.get('L', 1.0)

72 super().__init__(alpha=1.0, dtype=np.dtype('float'), **kwargs)

75 def _eval_explicit_part(self, u, t, f_expl):

76 f_expl[:] = -2.0 / self.eps**2 * u * (1.0 - u) * (1.0 - 2.0 * u) - 6.0 * self.dw * u * (1.0 - u)

77 return f_expl

79 def eval_f(self, u, t):

80 """

81 Routine to evaluate the right-hand side of the problem.

83 Parameters

84 ----------

85 u : dtype_u

86 Current values of the numerical solution.

87 t : float

88 Current time of the numerical solution is computed.

90 Returns

91 -------

92 f : dtype_f

93 The right-hand side of the problem.

94 """

96 f = self.dtype_f(self.init)

98 f.impl[:] = self._eval_Laplacian(u, f.impl)

100 if self.spectral:

101 f.impl = -self.K2 * u

103 if self.eps > 0:

104 tmp = self.fft.backward(u)

105 tmp[:] = self._eval_explicit_part(tmp, t, tmp)

106 f.expl[:] = self.fft.forward(tmp)

108 else:

110 if self.eps > 0:

111 f.expl[:] = self._eval_explicit_part(u, t, f.expl)

113 self.work_counters['rhs']()

114 return f

116 def u_exact(self, t, **kwargs):

117 r"""

118 Routine to compute the exact solution at time :math:t.

120 Parameters

121 ----------

122 t : float

123 Time of the exact solution.

125 Returns

126 -------

127 me : dtype_u

128 The exact solution.

129 """

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

132 me = self.dtype_u(self.init, val=0.0)

133 if self.init_type == 'circle':

134 r2 = (self.X[0] - 0.5) ** 2 + (self.X[1] - 0.5) ** 2

135 if self.spectral:

136 tmp = 0.5 * (1.0 + self.xp.tanh((self.radius - self.xp.sqrt(r2)) / (np.sqrt(2) * self.eps)))

137 me[:] = self.fft.forward(tmp)

138 else:

139 me[:] = 0.5 * (1.0 + self.xp.tanh((self.radius - self.xp.sqrt(r2)) / (np.sqrt(2) * self.eps)))

140 elif self.init_type == 'circle_rand':

141 ndim = len(me.shape)

142 L = int(self.L[0])

143 # get random radii for circles/spheres

144 self.xp.random.seed(1)

145 lbound = 3.0 * self.eps

146 ubound = 0.5 - self.eps

147 rand_radii = (ubound - lbound) * self.xp.random.random_sample(size=tuple([L] * ndim)) + lbound

148 # distribute circles/spheres

149 tmp = newDistArray(self.fft, False)

150 if ndim == 2:

151 for i in range(0, L):

152 for j in range(0, L):

154 r2 = (self.X[0] + i - L + 0.5) ** 2 + (self.X[1] + j - L + 0.5) ** 2

155 # add this blob, shifted by 1 to avoid issues with adding up negative contributions

156 tmp += self.xp.tanh((rand_radii[i, j] - self.xp.sqrt(r2)) / (np.sqrt(2) * self.eps)) + 1

157 else:

158 raise NotImplementedError

159 # normalize to [0,1]

160 tmp *= 0.5

161 assert self.xp.all(tmp <= 1.0)

162 if self.spectral:

163 me[:] = self.fft.forward(tmp)

164 else:

165 self.xp.copyto(me, tmp)

166 else:

167 raise NotImplementedError('type of initial value not implemented, got %s' % self.init_type)

169 return me

172class allencahn_imex_timeforcing(allencahn_imex):

173 r"""

174 Example implementing the :math:N-dimensional Allen-Cahn equation with periodic boundary conditions :math:u \in [0, 1]^2

175 using time-dependent forcing

177 .. math::

178 \frac{\partial u}{\partial t} = \Delta u - \frac{2}{\varepsilon^2} u (1 - u) (1 - 2u)

179 - 6 d_w u (1 - u)

181 on a spatial domain :math:[-\frac{L}{2}, \frac{L}{2}]^2, driving force :math:d_w, and :math:N=2,3. Different initial

182 conditions can be used, for example, circles of the form

184 .. math::

185 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{(x_i-0.5)^2 + (y_j-0.5)^2}}{\sqrt{2}\varepsilon}\right),

187 for :math:i, j=0,..,N-1, where :math:N is the number of spatial grid points. For time-stepping, the problem is treated

188 *semi-implicitly*, i.e., the linear part is solved with Fast-Fourier Transform (FFT) using mpi4py-fft [1]_ and the nonlinear part in the right-hand

189 side will be treated explicitly.

190 """

192 def eval_f(self, u, t):

193 """

194 Routine to evaluate the right-hand side of the problem.

196 Parameters

197 ----------

198 u : dtype_u

199 Current values of the numerical solution.

200 t : float

201 Current time of the numerical solution is computed.

203 Returns

204 -------

205 f : dtype_f

206 The right-hand side of the problem.

207 """

209 f = self.dtype_f(self.init)

211 f.impl[:] = self._eval_Laplacian(u, f.impl)

213 if self.spectral:

215 tmp = newDistArray(self.fft, False)

216 tmp[:] = self.fft.backward(u, tmp)

218 if self.eps > 0:

219 tmpf = -2.0 / self.eps**2 * tmp * (1.0 - tmp) * (1.0 - 2.0 * tmp)

220 else:

221 tmpf = self.dtype_f(self.init, val=0.0)

223 # build sum over RHS without driving force

224 Rt_local = float(self.xp.sum(self.fft.backward(f.impl) + tmpf))

225 if self.comm is not None:

226 Rt_global = self.comm.allreduce(sendobj=Rt_local, op=MPI.SUM)

227 else:

228 Rt_global = Rt_local

230 # build sum over driving force term

231 Ht_local = float(self.xp.sum(6.0 * tmp * (1.0 - tmp)))

232 if self.comm is not None:

233 Ht_global = self.comm.allreduce(sendobj=Ht_local, op=MPI.SUM)

234 else:

235 Ht_global = Ht_local

237 # add/substract time-dependent driving force

238 if Ht_global != 0.0:

239 dw = Rt_global / Ht_global

240 else:

241 dw = 0.0

243 tmpf -= 6.0 * dw * tmp * (1.0 - tmp)

244 f.expl[:] = self.fft.forward(tmpf)

246 else:

248 if self.eps > 0:

249 f.expl = -2.0 / self.eps**2 * u * (1.0 - u) * (1.0 - 2.0 * u)

251 # build sum over RHS without driving force

252 Rt_local = float(self.xp.sum(f.impl + f.expl))

253 if self.comm is not None:

254 Rt_global = self.comm.allreduce(sendobj=Rt_local, op=MPI.SUM)

255 else:

256 Rt_global = Rt_local

258 # build sum over driving force term

259 Ht_local = float(self.xp.sum(6.0 * u * (1.0 - u)))

260 if self.comm is not None:

261 Ht_global = self.comm.allreduce(sendobj=Ht_local, op=MPI.SUM)

262 else:

263 Ht_global = Ht_local

265 # add/substract time-dependent driving force

266 if Ht_global != 0.0:

267 dw = Rt_global / Ht_global

268 else:

269 dw = 0.0

271 f.expl -= 6.0 * dw * u * (1.0 - u)

273 self.work_counters['rhs']()

274 return f