Coverage for pySDC/implementations/problem_classes/RayleighBenard.py: 75%

1import numpy as np

2from mpi4py import MPI

4from pySDC.implementations.problem_classes.generic_spectral import GenericSpectralLinear

5from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh

6from pySDC.core.convergence_controller import ConvergenceController

7from pySDC.core.hooks import Hooks

8from pySDC.implementations.convergence_controller_classes.check_convergence import CheckConvergence

9from pySDC.core.problem import WorkCounter

12class RayleighBenard(GenericSpectralLinear):

13 """

14 Rayleigh-Benard Convection is a variation of incompressible Navier-Stokes.

16 The equations we solve are

18 u_x + v_z = 0

19 T_t - kappa (T_xx + T_zz) = -uT_x - vT_z

20 u_t - nu (u_xx + u_zz) + p_x = -uu_x - vu_z

21 v_t - nu (v_xx + v_zz) + p_z - T = -uv_x - vv_z

23 with u the horizontal velocity, v the vertical velocity (in z-direction), T the temperature, p the pressure, indices

24 denoting derivatives, kappa=(Rayleigh * Prandtl)**(-1/2) and nu = (Rayleigh / Prandtl)**(-1/2). Everything on the left

25 hand side, that is the viscous part, the pressure gradient and the buoyancy due to temperature are treated

26 implicitly, while the non-linear convection part on the right hand side is integrated explicitly.

28 The domain, vertical boundary conditions and pressure gauge are

30 Omega = [0, 8) x (-1, 1)

31 T(z=+1) = 0

32 T(z=-1) = 2

33 u(z=+-1) = v(z=+-1) = 0

34 integral over p = 0

36 The spectral discretization uses FFT horizontally, implying periodic BCs, and an ultraspherical method vertically to

37 facilitate the Dirichlet BCs.

39 Parameters:

40 Prandtl (float): Prandtl number

41 Rayleigh (float): Rayleigh number

42 nx (int): Horizontal resolution

43 nz (int): Vertical resolution

44 BCs (dict): Can specify boundary conditions here

45 dealiasing (float): Dealiasing factor for evaluating the non-linear part

46 comm (mpi4py.Intracomm): Space communicator

47 """

49 dtype_u = mesh

50 dtype_f = imex_mesh

52 def __init__(

53 self,

54 Prandtl=1,

55 Rayleigh=2e6,

56 nx=256,

57 nz=64,

58 BCs=None,

59 dealiasing=3 / 2,

60 comm=None,

61 Lx=4,

62 Lz=1,

63 z0=0,

64 **kwargs,

65 ):

66 """

67 Constructor. `kwargs` are forwarded to parent class constructor.

69 Args:

70 Prandtl (float): Prandtl number

71 Rayleigh (float): Rayleigh number

72 nx (int): Resolution in x-direction

73 nz (int): Resolution in z direction

74 BCs (dict): Vertical boundary conditions

75 dealiasing (float): Dealiasing for evaluating the non-linear part in real space

76 comm (mpi4py.Intracomm): Space communicator

77 Lx (float): Horizontal length of the domain

78 Lz (float): Vertical length of the domain

79 z0 (float): Position of lower boundary

80 """

81 BCs = {} if BCs is None else BCs

82 BCs = {

83 'T_top': 0,

84 'T_bottom': 1,

85 'v_top': 0,

86 'v_bottom': 0,

87 'u_top': 0,

88 'u_bottom': 0,

89 'p_integral': 0,

90 **BCs,

91 }

92 if comm is None:

93 try:

94 from mpi4py import MPI

96 comm = MPI.COMM_WORLD

97 except ModuleNotFoundError:

98 pass

99 self._makeAttributeAndRegister(

100 'Prandtl',

101 'Rayleigh',

102 'nx',

103 'nz',

104 'BCs',

105 'dealiasing',

106 'comm',

107 'Lx',

108 'Lz',

109 'z0',

110 localVars=locals(),

111 readOnly=True,

112 )

113

114 bases = [

115 {'base': 'fft', 'N': nx, 'x0': 0, 'x1': self.Lx},

116 {'base': 'ultraspherical', 'N': nz, 'x0': self.z0, 'x1': self.Lz},

117 ]

118 components = ['u', 'v', 'T', 'p']

119 super().__init__(bases, components, comm=comm, **kwargs)

120

121 self.X, self.Z = self.get_grid()

122 self.Kx, self.Kz = self.get_wavenumbers()

123

124 # construct 2D matrices

125 Dzz = self.get_differentiation_matrix(axes=(1,), p=2)

126 Dz = self.get_differentiation_matrix(axes=(1,))

127 Dx = self.get_differentiation_matrix(axes=(0,))

128 Dxx = self.get_differentiation_matrix(axes=(0,), p=2)

129 Id = self.get_Id()

130

131 S1 = self.get_basis_change_matrix(p_out=0, p_in=1)

132 S2 = self.get_basis_change_matrix(p_out=0, p_in=2)

133

134 U01 = self.get_basis_change_matrix(p_in=0, p_out=1)

135 U12 = self.get_basis_change_matrix(p_in=1, p_out=2)

136 U02 = self.get_basis_change_matrix(p_in=0, p_out=2)

137

138 self.Dx = Dx

139 self.Dxx = Dxx

140 self.Dz = S1 @ Dz

141 self.Dzz = S2 @ Dzz

142

143 # compute rescaled Rayleigh number to extract viscosity and thermal diffusivity

144 Ra = Rayleigh / (max([abs(BCs['T_top'] - BCs['T_bottom']), np.finfo(float).eps]) * self.axes[1].L ** 3)

145 self.kappa = (Ra * Prandtl) ** (-1 / 2.0)

146 self.nu = (Ra / Prandtl) ** (-1 / 2.0)

147

148 # construct operators

149 L_lhs = {

150 'p': {'u': U01 @ Dx, 'v': Dz}, # divergence free constraint

151 'u': {'p': U02 @ Dx, 'u': -self.nu * (U02 @ Dxx + Dzz)},

152 'v': {'p': U12 @ Dz, 'v': -self.nu * (U02 @ Dxx + Dzz), 'T': -U02 @ Id},

153 'T': {'T': -self.kappa * (U02 @ Dxx + Dzz)},

154 }

155 self.setup_L(L_lhs)

156

157 # mass matrix

158 M_lhs = {i: {i: U02 @ Id} for i in ['u', 'v', 'T']}

159 self.setup_M(M_lhs)

160

161 # Prepare going from second (first for divergence free equation) derivative basis back to Chebychov-T

162 self.base_change = self._setup_operator({**{comp: {comp: S2} for comp in ['u', 'v', 'T']}, 'p': {'p': S1}})

163

164 # BCs

165 self.add_BC(

166 component='p', equation='p', axis=1, v=self.BCs['p_integral'], kind='integral', line=-1, scalar=True

167 )

168 self.add_BC(component='T', equation='T', axis=1, x=-1, v=self.BCs['T_bottom'], kind='Dirichlet', line=-1)

169 self.add_BC(component='T', equation='T', axis=1, x=1, v=self.BCs['T_top'], kind='Dirichlet', line=-2)

170 self.add_BC(component='v', equation='v', axis=1, x=1, v=self.BCs['v_top'], kind='Dirichlet', line=-1)

171 self.add_BC(component='v', equation='v', axis=1, x=-1, v=self.BCs['v_bottom'], kind='Dirichlet', line=-2)

172 self.remove_BC(component='v', equation='v', axis=1, x=-1, kind='Dirichlet', line=-2, scalar=True)

173 self.add_BC(component='u', equation='u', axis=1, v=self.BCs['u_top'], x=1, kind='Dirichlet', line=-2)

174 self.add_BC(

175 component='u',

176 equation='u',

177 axis=1,

178 v=self.BCs['u_bottom'],

179 x=-1,

180 kind='Dirichlet',

181 line=-1,

182 )

183

184 # eliminate Nyquist mode if needed

185 if nx % 2 == 0:

186 Nyquist_mode_index = self.axes[0].get_Nyquist_mode_index()

187 for component in self.components:

188 self.add_BC(

189 component=component, equation=component, axis=0, kind='Nyquist', line=int(Nyquist_mode_index), v=0

190 )

191 self.setup_BCs()

192

193 self.work_counters['rhs'] = WorkCounter()

194

195 def eval_f(self, u, *args, **kwargs):

196 f = self.f_init

197

198 if self.spectral_space:

199 u_hat = u.copy()

200 else:

201 u_hat = self.transform(u)

202

203 f_impl_hat = self.u_init_forward

204

205 Dz = self.Dz

206 Dx = self.Dx

207

208 iu, iv, iT, ip = self.index(['u', 'v', 'T', 'p'])

209

210 # evaluate implicit terms

211 if not hasattr(self, '_L_T_base'):

212 self._L_T_base = self.base_change @ self.L

213 f_impl_hat = -(self._L_T_base @ u_hat.flatten()).reshape(u_hat.shape)

214

215 if self.spectral_space:

216 f.impl[:] = f_impl_hat

217 else:

218 f.impl[:] = self.itransform(f_impl_hat).real

219

220 # -------------------------------------------

221 # treat convection explicitly with dealiasing

222

223 # start by computing derivatives

224 if not hasattr(self, '_Dx_expanded') or not hasattr(self, '_Dz_expanded'):

225 self._Dx_expanded = self._setup_operator({'u': {'u': Dx}, 'v': {'v': Dx}, 'T': {'T': Dx}, 'p': {}})

226 self._Dz_expanded = self._setup_operator({'u': {'u': Dz}, 'v': {'v': Dz}, 'T': {'T': Dz}, 'p': {}})

227 Dx_u_hat = (self._Dx_expanded @ u_hat.flatten()).reshape(u_hat.shape)

228 Dz_u_hat = (self._Dz_expanded @ u_hat.flatten()).reshape(u_hat.shape)

229

230 padding = (self.dealiasing, self.dealiasing)

231 Dx_u_pad = self.itransform(Dx_u_hat, padding=padding).real

232 Dz_u_pad = self.itransform(Dz_u_hat, padding=padding).real

233 u_pad = self.itransform(u_hat, padding=padding).real

234

235 fexpl_pad = self.xp.zeros_like(u_pad)

236 fexpl_pad[iu][:] = -(u_pad[iu] * Dx_u_pad[iu] + u_pad[iv] * Dz_u_pad[iu])

237 fexpl_pad[iv][:] = -(u_pad[iu] * Dx_u_pad[iv] + u_pad[iv] * Dz_u_pad[iv])

238 fexpl_pad[iT][:] = -(u_pad[iu] * Dx_u_pad[iT] + u_pad[iv] * Dz_u_pad[iT])

239

240 if self.spectral_space:

241 f.expl[:] = self.transform(fexpl_pad, padding=padding)

242 else:

243 f.expl[:] = self.itransform(self.transform(fexpl_pad, padding=padding)).real

244

245 self.work_counters['rhs']()

246 return f

247

248 def u_exact(self, t=0, noise_level=1e-3, seed=99):

249 assert t == 0

250 assert (

251 self.BCs['v_top'] == self.BCs['v_bottom']

252 ), 'Initial conditions are only implemented for zero velocity gradient'

253

254 me = self.spectral.u_init

255 iu, iv, iT, ip = self.index(['u', 'v', 'T', 'p'])

256

257 # linear temperature gradient

258 for comp in ['T', 'v', 'u']:

259 a = (self.BCs[f'{comp}_top'] - self.BCs[f'{comp}_bottom']) / self.Lz

260 b = self.BCs[f'{comp}_bottom'] - a * self.z0

261 me[self.index(comp)] = a * self.Z + b

262

263 # perturb slightly

264 rng = self.xp.random.default_rng(seed=seed)

265

266 noise = self.spectral.u_init

267 noise[iT] = rng.random(size=me[iT].shape)

268

269 me[iT] += noise[iT].real * noise_level * (self.Z - self.z0) * (self.Z - self.z0 + self.Lz)

270

271 if self.spectral_space:

272 me_hat = self.spectral.u_init_forward

273 me_hat[:] = self.transform(me)

274 return me_hat

275 else:

276 return me

277

278 def apply_BCs(self, sol):

279 """

280 Enforce the Dirichlet BCs at the top and bottom for arbitrary solution.

281 The function modifies the last two modes of u, v, and T in order to achieve this.

282 Note that the pressure is not modified here and the Nyquist mode is not altered either.

283

284 Args:

285 sol: Some solution that does not need to enforce boundary conditions

286

287 Returns:

288 Modified version of the solution that satisfies Dirichlet BCs.

289 """

290 ultraspherical = self.spectral.axes[-1]

291

292 if self.spectral_space:

293 sol_half_hat = self.itransform(sol, axes=(-2,))

294 else:

295 sol_half_hat = self.transform(sol, axes=(-1,))

296

297 BC_bottom = ultraspherical.get_BC(x=-1, kind='dirichlet')

298 BC_top = ultraspherical.get_BC(x=1, kind='dirichlet')

299

300 M = np.array([BC_top[-2:], BC_bottom[-2:]])

301 M_I = np.linalg.inv(M)

302 rhs = np.empty((2, self.nx), dtype=complex)

303 for component in ['u', 'v', 'T']:

304 i = self.index(component)

305 rhs[0] = self.BCs[f'{component}_top'] - self.xp.sum(sol_half_hat[i, :, :-2] * BC_top[:-2], axis=1)

306 rhs[1] = self.BCs[f'{component}_bottom'] - self.xp.sum(sol_half_hat[i, :, :-2] * BC_bottom[:-2], axis=1)

307

308 BC_vals = M_I @ rhs

309

310 sol_half_hat[i, :, -2:] = BC_vals.T

311

312 if self.spectral_space:

313 return self.transform(sol_half_hat, axes=(-2,))

314 else:

315 return self.itransform(sol_half_hat, axes=(-1,))

316

317 def get_fig(self): # pragma: no cover

318 """

319 Get a figure suitable to plot the solution of this problem

320

321 Returns

322 -------

323 self.fig : matplotlib.pyplot.figure.Figure

324 """

325 import matplotlib.pyplot as plt

326 from mpl_toolkits.axes_grid1 import make_axes_locatable

327

328 plt.rcParams['figure.constrained_layout.use'] = True

329 self.fig, axs = plt.subplots(2, 1, sharex=True, sharey=True, figsize=((10, 5)))

330 self.cax = []

331 divider = make_axes_locatable(axs[0])

332 self.cax += [divider.append_axes('right', size='3%', pad=0.03)]

333 divider2 = make_axes_locatable(axs[1])

334 self.cax += [divider2.append_axes('right', size='3%', pad=0.03)]

335 return self.fig

336

337 def plot(self, u, t=None, fig=None, quantity='T'): # pragma: no cover

338 r"""

339 Plot the solution.

340

341 Parameters

342 ----------

343 u : dtype_u

344 Solution to be plotted

345 t : float

346 Time to display at the top of the figure

347 fig : matplotlib.pyplot.figure.Figure

348 Figure with the same structure as a figure generated by `self.get_fig`. If none is supplied, a new figure will be generated.

349 quantity : (str)

350 quantity you want to plot

351

352 Returns

353 -------

354 None

355 """

356 fig = self.get_fig() if fig is None else fig

357 axs = fig.axes

358

359 imV = axs[1].pcolormesh(self.X, self.Z, self.compute_vorticity(u).real)

360

361 if self.spectral_space:

362 u = self.itransform(u)

363

364 imT = axs[0].pcolormesh(self.X, self.Z, u[self.index(quantity)].real)

365

366 for i, label in zip([0, 1], [rf'${quantity}$', 'vorticity']):

367 axs[i].set_aspect(1)

368 axs[i].set_title(label)

369

370 if t is not None:

371 fig.suptitle(f't = {t:.2f}')

372 axs[1].set_xlabel(r'$x$')

373 axs[1].set_ylabel(r'$z$')

374 fig.colorbar(imT, self.cax[0])

375 fig.colorbar(imV, self.cax[1])

376

377 def compute_vorticity(self, u):

378 if self.spectral_space:

379 u_hat = u.copy()

380 else:

381 u_hat = self.transform(u)

382

383 Dz = self.Dz

384 Dx = self.Dx

385 iu, iv = self.index(['u', 'v'])

386

387 vorticity_hat = self.spectral.u_init_forward

388 vorticity_hat[0] = (Dx * u_hat[iv].flatten() + Dz @ u_hat[iu].flatten()).reshape(u_hat[iu].shape)

389 return self.itransform(vorticity_hat)[0].real

390

391 def getOutputFile(self, fileName):

392 from pySDC.helpers.fieldsIO import Rectilinear

393

394 self.setUpFieldsIO()

395

396 coords = [me.get_1dgrid() for me in self.spectral.axes]

397 assert np.allclose([len(me) for me in coords], self.spectral.global_shape[1:])

398

399 fOut = Rectilinear(np.float64, fileName=fileName)

400 fOut.setHeader(nVar=len(self.components) + 1, coords=coords)

401 fOut.initialize()

402 return fOut

403

404 def processSolutionForOutput(self, u):

405 vorticity = self.compute_vorticity(u)

406

407 if self.spectral_space:

408 u_real = self.itransform(u).real

409 else:

410 u_real = u.real

411

412 me = np.empty(shape=(u_real.shape[0] + 1, *vorticity.shape))

413 me[:-1] = u_real

414 me[-1] = vorticity

415 return me

416

417 def compute_Nusselt_numbers(self, u):

418 """

419 Compute the various versions of the Nusselt number. This reflects the type of heat transport.

420 If the Nusselt number is equal to one, it indicates heat transport due to conduction. If it is larger,

421 advection is present.

422 Computing the Nusselt number at various places can be used to check the code.

423

424 Args:

425 u: The solution you want to compute the Nusselt numbers of

426

427 Returns:

428 dict: Nusselt number averaged over the entire volume and horizontally averaged at the top and bottom.

429 """

430 iv, iT = self.index(['v', 'T'])

431 zAxis = self.spectral.axes[-1]

432

433 if self.spectral_space:

434 u_hat = u.copy()

435 else:

436 u_hat = self.transform(u)

437

438 DzT_hat = (self.Dz @ u_hat[iT].flatten()).reshape(u_hat[iT].shape)

439

440 # compute vT with dealiasing

441 padding = (self.dealiasing, self.dealiasing)

442 u_pad = self.itransform(u_hat, padding=padding).real

443 _me = self.xp.zeros_like(u_pad)

444 _me[0] = u_pad[iv] * u_pad[iT]

445 vT_hat = self.transform(_me, padding=padding)[0]

446

447 if not hasattr(self, '_zInt'):

448 self._zInt = zAxis.get_integration_matrix()

449

450 nusselt_hat = (vT_hat / self.kappa - DzT_hat) * self.axes[-1].L

451

452 # get coefficients for evaluation on the boundary

453 top = zAxis.get_BC(kind='Dirichlet', x=1)

454 bot = zAxis.get_BC(kind='Dirichlet', x=-1)

455

456 integral_V = 0

457 if self.comm.rank == 0:

458

459 integral_z = (self._zInt @ nusselt_hat[0]).real

460 integral_z[0] = zAxis.get_integration_constant(integral_z, axis=-1)

461 integral_V = ((top - bot) * integral_z).sum() * self.axes[0].L / self.nx

462

463 Nusselt_V = self.comm.bcast(integral_V / self.spectral.V, root=0)

464 Nusselt_t = self.comm.bcast(self.xp.sum(nusselt_hat.real[0] * top, axis=-1) / self.nx, root=0)

465 Nusselt_b = self.comm.bcast(self.xp.sum(nusselt_hat.real[0] * bot, axis=-1) / self.nx, root=0)

466

467 return {

468 'V': Nusselt_V,

469 't': Nusselt_t,

470 'b': Nusselt_b,

471 }

472

473 def compute_viscous_dissipation(self, u):

474 iu, iv = self.index(['u', 'v'])

475

476 Lap_u_hat = self.spectral.u_init_forward

477

478 if self.spectral_space:

479 u_hat = u.copy()

480 else:

481 u_hat = self.transform(u)

482 Lap_u_hat[iu] = ((self.Dzz + self.Dxx) @ u_hat[iu].flatten()).reshape(u_hat[iu].shape)

483 Lap_u_hat[iv] = ((self.Dzz + self.Dxx) @ u_hat[iv].flatten()).reshape(u_hat[iu].shape)

484 Lap_u = self.itransform(Lap_u_hat)

485

486 return abs(u[iu] * Lap_u[iu] + u[iv] * Lap_u[iv])

487

488 def compute_buoyancy_generation(self, u):

489 if self.spectral_space:

490 u = self.itransform(u)

491 iv, iT = self.index(['v', 'T'])

492 return abs(u[iv] * self.Rayleigh * u[iT])

493

494

495class CFLLimit(ConvergenceController):

496

497 def dependencies(self, controller, *args, **kwargs):

498 from pySDC.implementations.hooks.log_step_size import LogStepSize

499

500 controller.add_hook(LogCFL)

501 controller.add_hook(LogStepSize)

502

503 def setup_status_variables(self, controller, **kwargs):

504 """

505 Add the embedded error variable to the error function.

506

507 Args:

508 controller (pySDC.Controller): The controller

509 """

510 self.add_status_variable_to_level('CFL_limit')

511

512 def setup(self, controller, params, description, **kwargs):

513 """

514 Define default parameters here.

515

516 Default parameters are:

517 - control_order (int): The order relative to other convergence controllers

518 - dt_max (float): maximal step size

519 - dt_min (float): minimal step size

520

521 Args:

522 controller (pySDC.Controller): The controller

523 params (dict): The params passed for this specific convergence controller

524 description (dict): The description object used to instantiate the controller

525

526 Returns:

527 (dict): The updated params dictionary

528 """

529 defaults = {

530 "control_order": -50,

531 "dt_max": np.inf,

532 "dt_min": 0,

533 "cfl": 0.4,

534 }

535 return {**defaults, **super().setup(controller, params, description, **kwargs)}

536

537 @staticmethod

538 def compute_max_step_size(P, u):

539 grid_spacing_x = P.X[1, 0] - P.X[0, 0]

540

541 cell_wallz = P.xp.zeros(P.nz + 1)

542 cell_wallz[0] = P.Lz

543 cell_wallz[-1] = 0

544 cell_wallz[1:-1] = (P.Z[0, :-1] + P.Z[0, 1:]) / 2

545 grid_spacing_z = cell_wallz[:-1] - cell_wallz[1:]

546

547 iu, iv = P.index(['u', 'v'])

548

549 if P.spectral_space:

550 u = P.itransform(u)

551

552 max_step_size_x = P.xp.min(grid_spacing_x / P.xp.abs(u[iu]))

553 max_step_size_z = P.xp.min(grid_spacing_z / P.xp.abs(u[iv]))

554 max_step_size = min([max_step_size_x, max_step_size_z])

555

556 if hasattr(P, 'comm'):

557 max_step_size = P.comm.allreduce(max_step_size, op=MPI.MIN)

558 return float(max_step_size)

559

560 def get_new_step_size(self, controller, step, **kwargs):

561 if not CheckConvergence.check_convergence(step):

562 return None

563

564 L = step.levels[0]

565 P = step.levels[0].prob

566

567 L.sweep.compute_end_point()

568 max_step_size = self.compute_max_step_size(P, L.uend)

569

570 L.status.CFL_limit = self.params.cfl * max_step_size

571

572 dt_new = L.status.dt_new if L.status.dt_new else max([self.params.dt_max, L.params.dt])

573 L.status.dt_new = min([dt_new, self.params.cfl * max_step_size])

574 L.status.dt_new = max([self.params.dt_min, L.status.dt_new])

575

576 self.log(f'dt max: {max_step_size:.2e} -> New step size: {L.status.dt_new:.2e}', step)

577

578

579class LogCFL(Hooks):

580

581 def post_step(self, step, level_number):

582 """

583 Record CFL limit.

584

585 Args:

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

587 level_number (int): the current level number

588

589 Returns:

590 None

591 """

592 super().post_step(step, level_number)

593

594 L = step.levels[level_number]

595

596 self.add_to_stats(

597 process=step.status.slot,

598 time=L.time + L.dt,

599 level=L.level_index,

600 iter=step.status.iter,

601 sweep=L.status.sweep,

602 type='CFL_limit',

603 value=L.status.CFL_limit,

604 )

605

606

607class LogAnalysisVariables(Hooks):

608

609 def post_step(self, step, level_number):

610 """

611 Record Nusselt numbers.

612

613 Args:

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

615 level_number (int): the current level number

616

617 Returns:

618 None

619 """

620 super().post_step(step, level_number)

621

622 L = step.levels[level_number]

623 P = L.prob

624

625 L.sweep.compute_end_point()

626 Nusselt = P.compute_Nusselt_numbers(L.uend)

627 buoyancy_production = P.compute_buoyancy_generation(L.uend)

628 viscous_dissipation = P.compute_viscous_dissipation(L.uend)

629

630 for key, value in zip(

631 ['Nusselt', 'buoyancy_production', 'viscous_dissipation'],

632 [Nusselt, buoyancy_production, viscous_dissipation],

633 ):

634 self.add_to_stats(

635 process=step.status.slot,

636 time=L.time + L.dt,

637 level=L.level_index,

638 iter=step.status.iter,

639 sweep=L.status.sweep,

640 type=key,

641 value=value,

642 )