[pyfr] 72/88: DTS cleanups and convergence monitoring.
Ghislain Vaillant
ghisvail-guest at moszumanska.debian.org
Wed Nov 16 12:05:31 UTC 2016
This is an automated email from the git hooks/post-receive script.
ghisvail-guest pushed a commit to branch master
in repository pyfr.
commit 5d1b312cf2b6776d10bb9d22bced8809004bb8e1
Author: Niki Loppi <nal15 at imperial.ac.uk>
Date: Fri Aug 19 10:56:16 2016 +0100
DTS cleanups and convergence monitoring.
---
pyfr/integrators/dual/base.py | 6 +--
pyfr/integrators/dual/controllers.py | 50 +++++++++++++++--
pyfr/integrators/dual/pseudosteppers.py | 96 ++++++++++++++++-----------------
pyfr/integrators/dual/steppers.py | 10 ----
4 files changed, 96 insertions(+), 66 deletions(-)
diff --git a/pyfr/integrators/dual/base.py b/pyfr/integrators/dual/base.py
index dc5f266..b9de5a0 100644
--- a/pyfr/integrators/dual/base.py
+++ b/pyfr/integrators/dual/base.py
@@ -11,7 +11,7 @@ class BaseDualIntegrator(BaseIntegrator):
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
- self._dtau = self.cfg.getfloat('solver-time-integrator', 'dtau')
+ self._dtau = self.cfg.getfloat('solver-time-integrator', 'pseudo-dt')
self.dtaumin = 1.0e-12
@property
@@ -29,7 +29,3 @@ class BaseDualIntegrator(BaseIntegrator):
@abstractmethod
def finalise_step(self, currsoln):
pass
-
- @abstractmethod
- def _point_implicit_coeff(self, dt, stepper_coeff):
- pass
diff --git a/pyfr/integrators/dual/controllers.py b/pyfr/integrators/dual/controllers.py
index cd6a9f7..ea21e0f 100644
--- a/pyfr/integrators/dual/controllers.py
+++ b/pyfr/integrators/dual/controllers.py
@@ -1,6 +1,10 @@
# -*- coding: utf-8 -*-
+import math
+
from pyfr.integrators.dual.base import BaseDualIntegrator
+from pyfr.mpiutil import get_comm_rank_root, get_mpi
+from pyfr.util import memoize
class BaseDualController(BaseDualIntegrator):
@@ -41,22 +45,62 @@ class DualNoneController(BaseDualController):
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
- self._niters = self.cfg.getint('solver-time-integrator', 'niters')
+ sect = 'solver-time-integrator'
+
+ self._maxniters = self.cfg.getint(sect, 'pseudo-niters-max')
+ self._minniters = self.cfg.getint(sect, 'pseudo-niters-min')
+
+ if self._maxniters < self._minniters:
+ raise ValueError('The maximum number of pseudo-iterations must '
+ 'be greater than or equal to the minimum')
+
+ self._pseudo_aresid = self.cfg.getfloat(sect, 'pseudo-aresid')
+ self._pseudo_rresid = self.cfg.getfloat(sect, 'pseudo-rresid')
def advance_to(self, t):
if t < self.tcurr:
raise ValueError('Advance time is in the past')
while self.tcurr < t:
- for i in range(self._niters):
+ for i in range(self._maxniters):
dt = max(min(t - self.tcurr, self._dt), self.dtmin)
dtau = max(min(t - self.tcurr, self._dtau), self.dtaumin)
# Take the step
- self._idxcurr = self.step(self.tcurr, dt, dtau)
+ self._idxcurr, self._idxprev = self.step(self.tcurr, dt, dtau)
+
+ # Activate convergence monitoring after pseudo-niters-min
+ if i >= self._minniters - 1:
+ # Subtract the current solution from the previous solution
+ self._add(-1.0/dtau, self._idxprev, 1.0/dtau, self._idxcurr)
+
+ # Compute the normalised residual and check for convergence
+ if self._resid(self._idxprev, self._idxcurr) < 1.0:
+ break
# Update the dual-time stepping banks (n+1 => n, n => n-1)
self.finalise_step(self._idxcurr)
# We are not adaptive, so accept every step
self._accept_step(dt, self._idxcurr)
+
+ @memoize
+ def _get_errest_kerns(self):
+ return self._get_kernels('errest', nargs=3, norm='uniform')
+
+ def _resid(self, x, y):
+ comm, rank, root = get_comm_rank_root()
+
+ # Get an errest kern to compute the square of the maximum residual
+ errest = self._get_errest_kerns()
+
+ # Prepare and run the kernel
+ self._prepare_reg_banks(x, y, y)
+ self._queue % errest(self._pseudo_aresid, self._pseudo_rresid)
+
+ # Reduce locally (element types) and globally (MPI ranks)
+ rl = max(errest.retval)
+ rg = comm.allreduce(rl, op=get_mpi('max'))
+
+ # Normalise
+ return math.sqrt(rg)
diff --git a/pyfr/integrators/dual/pseudosteppers.py b/pyfr/integrators/dual/pseudosteppers.py
index 282dbb9..03861d5 100644
--- a/pyfr/integrators/dual/pseudosteppers.py
+++ b/pyfr/integrators/dual/pseudosteppers.py
@@ -10,19 +10,21 @@ class BaseDualPseudoStepper(BaseDualIntegrator):
stats.set('solver-time-integrator', 'nsteps', self.nsteps)
stats.set('solver-time-integrator', 'nfevals', self._stepper_nfevals)
- def _add_dual_source(self, dt, rhs, currsoln):
- # Coefficients for the dual-time source term
- coeffs = [c/dt for c in self._dual_time_source]
+ def _add_with_dts(self, *args, c):
+ vals, regs = list(args[::2]), list(args[1::2])
- # Get a suitable set of axnpby kernels
- axnpby = self._get_axnpby_kerns(len(coeffs) + 1,
- subdims=self._subdims)
+ # Coefficients for the dual-time source term
+ svals = [c*sc for sc in self._dual_time_source]
- # Prepare the matrix banks
- self._prepare_reg_banks(rhs, currsoln, *self._source_regidx)
+ # Normal addition
+ axnpby = self._get_axnpby_kerns(len(vals))
+ self._prepare_reg_banks(*regs)
+ self._queue % axnpby(*vals)
- # Bind and run the axnpby kernels
- self._queue % axnpby(1.0, *coeffs)
+ # Source addition
+ axnpby = self._get_axnpby_kerns(len(svals) + 1, subdims=self._subdims)
+ self._prepare_reg_banks(regs[0], self._idxcurr, *self._source_regidx)
+ self._queue % axnpby(1, *svals)
def finalise_step(self, currsoln):
add = self._add
@@ -52,15 +54,18 @@ class DualPseudoEulerStepper(BaseDualPseudoStepper):
return 1
def step(self, t, dt, dtau):
- add, add_dual_source = self._add, self._add_dual_source
+ add, add_with_dts = self._add, self._add_with_dts
rhs = self.system.rhs
- ut, f = self._stepper_regidx
+ r0, r1 = self._stepper_regidx
+ rat = dtau / dt
+
+ if r0 != self._idxcurr:
+ r0, r1 = r1, r0
- rhs(t, ut, f)
- add_dual_source(dt, f, ut)
- add(1.0, ut, dtau, f)
+ rhs(t, r0, r1)
+ add_with_dts(0, r1, 1, r0, dtau, r1, c=rat)
- return ut
+ return r1, r0
class DualPseudoTVDRK3Stepper(BaseDualPseudoStepper):
@@ -79,8 +84,9 @@ class DualPseudoTVDRK3Stepper(BaseDualPseudoStepper):
return 3
def step(self, t, dt, dtau):
- add, add_dual_source = self._add, self._add_dual_source
+ add, add_with_dts = self._add, self._add_with_dts
rhs = self.system.rhs
+ rat = dtau / dt
# Get the bank indices for pseudo-registers (n+1,m; n+1,m+1; rhs),
# where m = pseudo-time and n = real-time
@@ -91,25 +97,22 @@ class DualPseudoTVDRK3Stepper(BaseDualPseudoStepper):
r0, r1 = r1, r0
# First stage;
- # r2 = -∇·f(r0) - dQ/dt; r1 = r0 + dtau*r2
+ # r2 = -∇·f(r0); r1 = r0 + dtau*r2 - dtau*dQ/dt;
rhs(t, r0, r2)
- add_dual_source(dt, r2, r0)
- add(0.0, r1, 1.0, r0, dtau, r2)
+ add_with_dts(0, r1, 1, r0, dtau, r2, c=rat)
# Second stage;
- # r2 = -∇·f(r1) - dQ/dt; r1 = 0.75*r0 + 0.25*r1 + 0.25*dtau*r2
+ # r2 = -∇·f(r1); r1 = 3/4*r0 + 1/4*r1 + 1/4*dtau*r2 - dtau/4*dQ/dt
rhs(t, r1, r2)
- add_dual_source(dt, r2, r0)
- add(0.25, r1, 0.75, r0, dtau/4.0, r2)
+ add_with_dts(1/4, r1, 3/4, r0, dtau/4, r2, c=rat/4)
# Third stage;
- # r2 = -∇·f(r1) - dQ/dt; r1 = 1.0/3.0*r0 + 2.0/3.0*r1 + 2.0/3.0*dtau*r2
+ # r2 = -∇·f(r1); r1 = 1/3*r0 + 2/3*r1 + 2/3*dtau*r2 - 2/3*dtau*dQ/dt
rhs(t, r1, r2)
- add_dual_source(dt, r2, r1)
- add(2.0/3.0, r1, 1.0/3.0, r0, 2.0*dtau/3.0, r2)
+ add_with_dts(2/3, r1, 1/3, r0, 2*dtau/3, r2, c=2*rat/3)
# Return the index of the bank containing u(n+1,m+1)
- return r1
+ return r1, r0
class DualPseudoRK4Stepper(BaseDualPseudoStepper):
@@ -128,8 +131,9 @@ class DualPseudoRK4Stepper(BaseDualPseudoStepper):
return 4
def step(self, t, dt, dtau):
- add, add_dual_source = self._add, self._add_dual_source
- rhs, pi_coeff = self.system.rhs, self._point_implicit_coeff
+ add, add_with_dts = self._add, self._add_with_dts
+ rhs = self.system.rhs
+ rat = dtau / dt
# Get the bank indices for pseudo-registers (n+1,m; n+1,m+1; rhs),
# where m = pseudo-time and n = real-time
@@ -139,39 +143,35 @@ class DualPseudoRK4Stepper(BaseDualPseudoStepper):
if r0 != self._idxcurr:
r0, r1 = r1, r0
- # First stage; r1 = -∇·f(r0) - dQ/dt
+ # First stage; r1 = -∇·f(r0)
rhs(t, r0, r1)
- add_dual_source(dt, r1, r0)
- # Second stage; r2 = r0 + dtau/2*r1; r2 = -∇·f(r2) - dQ/dt
- add(0.0, r2, 1.0, r0, dtau/2.0, r1)
+ # Second stage; r2 = r0 + dtau/2*r1 - dtau/2*dQ/dt; r2 = -∇·f(r2)
+ add_with_dts(0, r2, 1, r0, dtau/2, r1, c=rat/2)
rhs(t, r2, r2)
- add_dual_source(dt, r2, r0)
# As no subsequent stages depend on the first stage we can
# reuse its register to start accumulating the solution with
- # r1 = r0 + pi_coeff*r1 + pi_coeff*r2
- add(pi_coeff(dt, dtau/6.0), r1, 1.0, r0, pi_coeff(dt, dtau/3.0), r2)
+ # r1 = r0 + dtau/6*r1 + dtau/3*r2 -(dtau/6+dtau/3)*dQ/dt
+ add_with_dts(dtau/6, r1, 1, r0, dtau/3, r2, c=(1/6+1/3)*rat)
# Third stage; here we reuse the r2 register
- # r2 = r0 + dtau/2*r2
+ # r2 = r0 + dtau/2*r2 - dtau/2*dQ/dt
# r2 = -∇·f(r2)
- add(dtau/2.0, r2, 1.0, r0)
+ add_with_dts(dtau/2, r2, 1, r0, c=rat/2)
rhs(t, r2, r2)
- add_dual_source(dt, r2, r0)
- # Accumulate; r1 = r1 + pi_coeff*r2
- add(1.0, r1, pi_coeff(dt, dtau/3.0), r2)
+ # Accumulate; r1 = r1 + dtau/3*r2 - dtau/3*dQ/dt
+ add_with_dts(1, r1, dtau/3, r2, c=rat/3)
# Fourth stage; again we reuse r2
- # r2 = r0 + dtau*r2
- # r2 = -∇·f(r2) - dQ/dt
- add(dtau, r2, 1.0, r0)
+ # r2 = r0 + dtau*r2 - dtau*dQ/dt
+ # r2 = -∇·f(r2)
+ add_with_dts(dtau, r2, 1, r0, c=rat)
rhs(t, r2, r2)
- add_dual_source(dt, r2, r0)
- # Final accumulation r1 = r1 + pi_coeff*r2 = u(n+1,m+1)
- add(1.0, r1, pi_coeff(dt, dtau/6.0), r2)
+ # Final accumulation r1 = r1 + dtau/6*r2 - dtau/6*dQ/dt = u(n+1,m+1)
+ add_with_dts(1, r1, dtau/6, r2, c=rat/6)
# Return the index of the bank containing u(n+1,m+1)
- return r1
+ return r1, r0
diff --git a/pyfr/integrators/dual/steppers.py b/pyfr/integrators/dual/steppers.py
index a623958..f63cf1f 100644
--- a/pyfr/integrators/dual/steppers.py
+++ b/pyfr/integrators/dual/steppers.py
@@ -27,10 +27,6 @@ class BDF2DualStepper(BaseDualStepper):
def _dual_time_source(self):
return [-1.5, 2.0, -0.5]
- # RK4 uses point-implicit treatment to save one matrix bank
- def _point_implicit_coeff(self, dt, stepper_coeff):
- return stepper_coeff/(1.0 + 1.5*stepper_coeff/dt)
-
class BDF3DualStepper(BaseDualStepper):
stepper_name = 'bdf3'
@@ -43,9 +39,6 @@ class BDF3DualStepper(BaseDualStepper):
def _dual_time_source(self):
return [-11.0/6.0, 3.0, -1.5, 1.0/3.0]
- def _point_implicit_coeff(self, dt, stepper_coeff):
- return stepper_coeff/(1.0 + 11.0*stepper_coeff/(6.0*dt))
-
class BackwardEulerDualStepper(BaseDualStepper):
stepper_name = 'backward-euler'
@@ -57,6 +50,3 @@ class BackwardEulerDualStepper(BaseDualStepper):
@property
def _dual_time_source(self):
return [-1.0, 1.0]
-
- def _point_implicit_coeff(self, dt, stepper_coeff):
- return stepper_coeff/(1.0 + stepper_coeff/dt)
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/pyfr.git
More information about the debian-science-commits
mailing list