# -*- coding: utf-8 -*-
import math
import numpy as np
from pyfr.integrators.std.base import BaseStdIntegrator
from pyfr.mpiutil import get_comm_rank_root, mpi
class BaseStdController(BaseStdIntegrator):
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
# Solution filtering frequency
self._fnsteps = self.cfg.getint('soln-filter', 'nsteps', '0')
# Stats on the most recent step
self.stepinfo = []
# Fire off any event handlers if not restarting
if not self.isrestart:
for csh in self.completed_step_handlers:
csh(self)
def _accept_step(self, dt, idxcurr, err=None):
self.tcurr += dt
self.nacptsteps += 1
self.nacptchain += 1
self.stepinfo.append((dt, 'accept', err))
self._idxcurr = idxcurr
# Filter
if self._fnsteps and self.nacptsteps % self._fnsteps == 0:
self.system.filt(idxcurr)
# Invalidate the solution cache
self._curr_soln = None
# Invalidate the solution gradients cache
self._curr_grad_soln = None
# Fire off any event handlers
for csh in self.completed_step_handlers:
csh(self)
# Abort if plugins request it
self._check_abort()
# Clear the step info
self.stepinfo = []
def _reject_step(self, dt, idxold, err=None):
if dt <= self.dtmin:
raise RuntimeError('Minimum sized time step rejected')
self.nacptchain = 0
self.nrjctsteps += 1
self.stepinfo.append((dt, 'reject', err))
self._idxcurr = idxold
[docs]class StdNoneController(BaseStdController):
controller_name = 'none'
@property
def controller_needs_errest(self):
return False
[docs] def advance_to(self, t):
if t < self.tcurr:
raise ValueError('Advance time is in the past')
while self.tcurr < t:
# Decide on the time step
dt = max(min(t - self.tcurr, self._dt), self.dtmin)
# Take the step
idxcurr = self.step(self.tcurr, dt)
# We are not adaptive, so accept every step
self._accept_step(dt, idxcurr)
[docs]class StdPIController(BaseStdController):
controller_name = 'pi'
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
sect = 'solver-time-integrator'
# Maximum time step
self.dtmax = self.cfg.getfloat(sect, 'dt-max', 1e2)
# Error tolerances
self._atol = self.cfg.getfloat(sect, 'atol')
self._rtol = self.cfg.getfloat(sect, 'rtol')
# Error norm
self._norm = self.cfg.get(sect, 'errest-norm', 'l2')
if self._norm not in {'l2', 'uniform'}:
raise ValueError('Invalid error norm')
# PI control values
self._alpha = self.cfg.getfloat(sect, 'pi-alpha', 0.58)
self._beta = self.cfg.getfloat(sect, 'pi-beta', 0.42)
# Estimate of previous error
self._errprev = 1.0
# Step size adjustment factors
self._saffac = self.cfg.getfloat(sect, 'safety-fact', 0.8)
self._maxfac = self.cfg.getfloat(sect, 'max-fact', 2.5)
self._minfac = self.cfg.getfloat(sect, 'min-fact', 0.3)
if not self._minfac < 1 <= self._maxfac:
raise ValueError('Invalid max-fact, min-fact')
@property
def controller_needs_errest(self):
return True
[docs] def _errest(self, rcurr, rprev, rerr):
comm, rank, root = get_comm_rank_root()
# Get a set of kernels to estimate the integration error
ekerns = self._get_reduction_kerns(rcurr, rprev, rerr, method='errest',
norm=self._norm)
# Bind the dynamic arguments
for kern in ekerns:
kern.bind(self._atol, self._rtol)
# Run the kernels
self.backend.run_kernels(ekerns, wait=True)
# Pseudo L2 norm
if self._norm == 'l2':
# Reduce locally (element types + field variables)
err = np.array([sum(v for k in ekerns for v in k.retval)])
# Reduce globally (MPI ranks)
comm.Allreduce(mpi.IN_PLACE, err, op=mpi.SUM)
# Normalise
err = math.sqrt(float(err) / self._gndofs)
# Uniform norm
else:
# Reduce locally (element types + field variables)
err = np.array([max(v for k in ekerns for v in k.retval)])
# Reduce globally (MPI ranks)
comm.Allreduce(mpi.IN_PLACE, err, op=mpi.MAX)
# Normalise
err = math.sqrt(float(err))
return err if not math.isnan(err) else 100
[docs] def advance_to(self, t):
if t < self.tcurr:
raise ValueError('Advance time is in the past')
# Constants
maxf = self._maxfac
minf = self._minfac
saff = self._saffac
sord = self.stepper_order
expa = self._alpha / sord
expb = self._beta / sord
while self.tcurr < t:
# Decide on the time step
dt = max(min(t - self.tcurr, self._dt, self.dtmax), self.dtmin)
# Take the step
idxcurr, idxprev, idxerr = self.step(self.tcurr, dt)
# Estimate the error
err = self._errest(idxcurr, idxprev, idxerr)
# Determine time step adjustment factor
fac = err**-expa * self._errprev**expb
fac = min(maxf, max(minf, saff*fac))
# Compute the size of the next step
self._dt = fac*dt
# Decide if to accept or reject the step
if err < 1.0:
self._errprev = err
self._accept_step(dt, idxcurr, err=err)
else:
self._reject_step(dt, idxprev, err=err)