[pyfr] 03/88: Added dual-time stepping.
Ghislain Vaillant
ghisvail-guest at moszumanska.debian.org
Wed Nov 16 12:05:23 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 3789e616a42b5cd63e485dc315d931af63f9ad14
Author: Niki Loppi <nal15 at imperial.ac.uk>
Date: Thu Feb 25 19:38:31 2016 +0000
Added dual-time stepping.
---
pyfr/integrators/__init__.py | 2 +-
pyfr/integrators/base.py | 19 +++-
pyfr/integrators/dual/base.py | 23 +++++
pyfr/integrators/dual/controllers.py | 57 ++++++++++-
pyfr/integrators/dual/pseudosteppers.py | 172 +++++++++++++++++++++++++++++++-
pyfr/integrators/dual/steppers.py | 57 ++++++++++-
pyfr/integrators/std/steppers.py | 15 ---
pyfr/solvers/euler/elements.py | 2 +
setup.py | 2 +
9 files changed, 329 insertions(+), 20 deletions(-)
diff --git a/pyfr/integrators/__init__.py b/pyfr/integrators/__init__.py
index 689ae2c..7e7651d 100644
--- a/pyfr/integrators/__init__.py
+++ b/pyfr/integrators/__init__.py
@@ -26,7 +26,7 @@ def get_integrator(backend, systemcls, rallocs, mesh, initsoln, cfg):
cc = subclass_where(BaseDualController, controller_name=cn)
pc = subclass_where(BaseDualPseudoStepper, pseudo_stepper_name=pn)
- sn = subclass_where(BaseDualStepper, stepper_name=sn)
+ sc = subclass_where(BaseDualStepper, stepper_name=sn)
bases = [(cn, cc), (pn, pc), (sn, sc)]
else:
diff --git a/pyfr/integrators/base.py b/pyfr/integrators/base.py
index 5bc3ccc..7a5fbdb 100644
--- a/pyfr/integrators/base.py
+++ b/pyfr/integrators/base.py
@@ -9,7 +9,7 @@ import numpy as np
from pyfr.inifile import Inifile
from pyfr.mpiutil import get_comm_rank_root, get_mpi
from pyfr.plugins import get_plugin
-from pyfr.util import proxylist
+from pyfr.util import memoize, proxylist
class BaseIntegrator(object, metaclass=ABCMeta):
@@ -73,6 +73,9 @@ class BaseIntegrator(object, metaclass=ABCMeta):
# Solution cache
self._curr_soln = None
+ # Add kernel cache
+ self._axnpby_kerns = {}
+
# Event handlers for advance_to
self.completed_step_handlers = proxylist(self._get_plugins())
@@ -128,6 +131,20 @@ class BaseIntegrator(object, metaclass=ABCMeta):
for reg, ix in zip(self._regs, bidxes):
reg.active = ix
+ @memoize
+ def _get_axnpby_kerns(self, n, subdims=None):
+ return self._get_kernels('axnpby', nargs=n, subdims=subdims)
+
+ def _add(self, *args):
+ # Get a suitable set of axnpby kernels
+ axnpby = self._get_axnpby_kerns(len(args) // 2)
+
+ # Bank indices are in odd-numbered arguments
+ self._prepare_reg_banks(*args[1::2])
+
+ # Bind and run the axnpby kernels
+ self._queue % axnpby(*args[::2])
+
def call_plugin_dt(self, dt):
ta = self.tlist
tb = deque(np.arange(self.tcurr, self.tend, dt).tolist())
diff --git a/pyfr/integrators/dual/base.py b/pyfr/integrators/dual/base.py
index 6c1036f..dc5f266 100644
--- a/pyfr/integrators/dual/base.py
+++ b/pyfr/integrators/dual/base.py
@@ -1,5 +1,7 @@
# -*- coding: utf-8 -*-
+from abc import abstractmethod
+
from pyfr.integrators.base import BaseIntegrator
@@ -9,4 +11,25 @@ class BaseDualIntegrator(BaseIntegrator):
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
+ self._dtau = self.cfg.getfloat('solver-time-integrator', 'dtau')
+ self.dtaumin = 1.0e-12
+
+ @property
+ def _stepper_regidx(self):
+ return self._regidx[:self._pseudo_stepper_nregs]
+
+ @property
+ def _source_regidx(self):
+ return self._regidx[self._pseudo_stepper_nregs:]
+
+ @abstractmethod
+ def _dual_time_source(self):
+ pass
+
+ @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 19bd585..cd6a9f7 100644
--- a/pyfr/integrators/dual/controllers.py
+++ b/pyfr/integrators/dual/controllers.py
@@ -4,4 +4,59 @@ from pyfr.integrators.dual.base import BaseDualIntegrator
class BaseDualController(BaseDualIntegrator):
- pass
+ 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 = []
+
+ def _accept_step(self, dt, idxcurr):
+ self.tcurr += dt
+ self.nacptsteps += 1
+ self.nacptchain += 1
+
+ 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
+
+ # Fire off any event handlers
+ self.completed_step_handlers(self)
+
+ @property
+ def nsteps(self):
+ return self.nacptsteps + self.nrjctsteps
+
+
+class DualNoneController(BaseDualController):
+ controller_name = 'none'
+
+ def __init__(self, *args, **kwargs):
+ super().__init__(*args, **kwargs)
+
+ self._niters = self.cfg.getint('solver-time-integrator', 'niters')
+
+ 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):
+ 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)
+
+ # 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)
diff --git a/pyfr/integrators/dual/pseudosteppers.py b/pyfr/integrators/dual/pseudosteppers.py
index 92667fc..282dbb9 100644
--- a/pyfr/integrators/dual/pseudosteppers.py
+++ b/pyfr/integrators/dual/pseudosteppers.py
@@ -4,4 +4,174 @@ from pyfr.integrators.dual.base import BaseDualIntegrator
class BaseDualPseudoStepper(BaseDualIntegrator):
- pass
+ def collect_stats(self, stats):
+ super().collect_stats(stats)
+
+ 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]
+
+ # Get a suitable set of axnpby kernels
+ axnpby = self._get_axnpby_kerns(len(coeffs) + 1,
+ subdims=self._subdims)
+
+ # Prepare the matrix banks
+ self._prepare_reg_banks(rhs, currsoln, *self._source_regidx)
+
+ # Bind and run the axnpby kernels
+ self._queue % axnpby(1.0, *coeffs)
+
+ def finalise_step(self, currsoln):
+ add = self._add
+ pnreg = self._pseudo_stepper_nregs
+
+ # Rotate the source registers to the right by one
+ self._regidx[pnreg:] = (self._source_regidx[-1:] +
+ self._source_regidx[:-1])
+
+ # Copy the current soln into the first source register
+ add(0.0, self._regidx[pnreg], 1.0, currsoln)
+
+
+class DualPseudoEulerStepper(BaseDualPseudoStepper):
+ pseudo_stepper_name = 'euler'
+
+ @property
+ def _stepper_nfevals(self):
+ return self.nsteps
+
+ @property
+ def _pseudo_stepper_nregs(self):
+ return 2
+
+ @property
+ def _pseudo_stepper_order(self):
+ return 1
+
+ def step(self, t, dt, dtau):
+ add, add_dual_source = self._add, self._add_dual_source
+ rhs = self.system.rhs
+ ut, f = self._stepper_regidx
+
+ rhs(t, ut, f)
+ add_dual_source(dt, f, ut)
+ add(1.0, ut, dtau, f)
+
+ return ut
+
+
+class DualPseudoTVDRK3Stepper(BaseDualPseudoStepper):
+ pseudo_stepper_name = 'tvd-rk3'
+
+ @property
+ def _stepper_nfevals(self):
+ return 3*self.nsteps
+
+ @property
+ def _pseudo_stepper_nregs(self):
+ return 3
+
+ @property
+ def _pseudo_stepper_order(self):
+ return 3
+
+ def step(self, t, dt, dtau):
+ add, add_dual_source = self._add, self._add_dual_source
+ rhs = self.system.rhs
+
+ # Get the bank indices for pseudo-registers (n+1,m; n+1,m+1; rhs),
+ # where m = pseudo-time and n = real-time
+ r0, r1, r2 = self._stepper_regidx
+
+ # Ensure r0 references the bank containing u(n+1,m)
+ if r0 != self._idxcurr:
+ r0, r1 = r1, r0
+
+ # First stage;
+ # r2 = -∇·f(r0) - dQ/dt; r1 = r0 + dtau*r2
+ rhs(t, r0, r2)
+ add_dual_source(dt, r2, r0)
+ add(0.0, r1, 1.0, r0, dtau, r2)
+
+ # Second stage;
+ # r2 = -∇·f(r1) - dQ/dt; r1 = 0.75*r0 + 0.25*r1 + 0.25*dtau*r2
+ rhs(t, r1, r2)
+ add_dual_source(dt, r2, r0)
+ add(0.25, r1, 0.75, r0, dtau/4.0, r2)
+
+ # Third stage;
+ # r2 = -∇·f(r1) - dQ/dt; r1 = 1.0/3.0*r0 + 2.0/3.0*r1 + 2.0/3.0*dtau*r2
+ 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)
+
+ # Return the index of the bank containing u(n+1,m+1)
+ return r1
+
+
+class DualPseudoRK4Stepper(BaseDualPseudoStepper):
+ pseudo_stepper_name = 'rk4'
+
+ @property
+ def _stepper_nfevals(self):
+ return 4*self.nsteps
+
+ @property
+ def _pseudo_stepper_nregs(self):
+ return 3
+
+ @property
+ def _pseudo_stepper_order(self):
+ 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
+
+ # Get the bank indices for pseudo-registers (n+1,m; n+1,m+1; rhs),
+ # where m = pseudo-time and n = real-time
+ r0, r1, r2 = self._stepper_regidx
+
+ # Ensure r0 references the bank containing u(n+1,m)
+ if r0 != self._idxcurr:
+ r0, r1 = r1, r0
+
+ # First stage; r1 = -∇·f(r0) - dQ/dt
+ 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)
+ 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)
+
+ # Third stage; here we reuse the r2 register
+ # r2 = r0 + dtau/2*r2
+ # r2 = -∇·f(r2)
+ add(dtau/2.0, r2, 1.0, r0)
+ 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)
+
+ # Fourth stage; again we reuse r2
+ # r2 = r0 + dtau*r2
+ # r2 = -∇·f(r2) - dQ/dt
+ add(dtau, r2, 1.0, r0)
+ 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)
+
+ # Return the index of the bank containing u(n+1,m+1)
+ return r1
diff --git a/pyfr/integrators/dual/steppers.py b/pyfr/integrators/dual/steppers.py
index 03a0761..a623958 100644
--- a/pyfr/integrators/dual/steppers.py
+++ b/pyfr/integrators/dual/steppers.py
@@ -4,4 +4,59 @@ from pyfr.integrators.dual.base import BaseDualIntegrator
class BaseDualStepper(BaseDualIntegrator):
- pass
+ def __init__(self, *args, **kwargs):
+ super().__init__(*args, **kwargs)
+
+ elementscls = self.system.elementscls
+ self._subdims = [elementscls.convarmap[self.system.ndims].index(v)
+ for v in elementscls.dualcoeffs[self.system.ndims]]
+
+ @property
+ def _stepper_nregs(self):
+ return self._pseudo_stepper_nregs + len(self._dual_time_source) - 1
+
+
+class BDF2DualStepper(BaseDualStepper):
+ stepper_name = 'bdf2'
+
+ @property
+ def _stepper_order(self):
+ return 2
+
+ @property
+ 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'
+
+ @property
+ def _stepper_order(self):
+ return 3
+
+ @property
+ 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'
+
+ @property
+ def _stepper_order(self):
+ return 1
+
+ @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)
diff --git a/pyfr/integrators/std/steppers.py b/pyfr/integrators/std/steppers.py
index c28c6f2..319cac9 100644
--- a/pyfr/integrators/std/steppers.py
+++ b/pyfr/integrators/std/steppers.py
@@ -1,7 +1,6 @@
# -*- coding: utf-8 -*-
from pyfr.integrators.std.base import BaseStdIntegrator
-from pyfr.util import memoize
class BaseStdStepper(BaseStdIntegrator):
@@ -17,20 +16,6 @@ class BaseStdStepper(BaseStdIntegrator):
stats.set('solver-time-integrator', 'nsteps', self.nsteps)
stats.set('solver-time-integrator', 'nfevals', self._stepper_nfevals)
- @memoize
- def _get_axnpby_kerns(self, n):
- return self._get_kernels('axnpby', nargs=n)
-
- def _add(self, *args):
- # Get a suitable set of axnpby kernels
- axnpby = self._get_axnpby_kerns(len(args) // 2)
-
- # Bank indices are in odd-numbered arguments
- self._prepare_reg_banks(*args[1::2])
-
- # Bind and run the axnpby kernels
- self._queue % axnpby(*args[::2])
-
class StdEulerStepper(BaseStdStepper):
stepper_name = 'euler'
diff --git a/pyfr/solvers/euler/elements.py b/pyfr/solvers/euler/elements.py
index fd94b1b..151f27e 100644
--- a/pyfr/solvers/euler/elements.py
+++ b/pyfr/solvers/euler/elements.py
@@ -12,6 +12,8 @@ class BaseFluidElements(object):
convarmap = {2: ['rho', 'rhou', 'rhov', 'E'],
3: ['rho', 'rhou', 'rhov', 'rhow', 'E']}
+ dualcoeffs = convarmap
+
visvarmap = {
2: {'density': ['rho'],
'velocity': ['u', 'v'],
diff --git a/setup.py b/setup.py
index fd7b355..3a89997 100755
--- a/setup.py
+++ b/setup.py
@@ -31,6 +31,8 @@ modules = [
'pyfr.backends.openmp',
'pyfr.backends.openmp.kernels',
'pyfr.integrators',
+ 'pyfr.integrators.dual',
+ 'pyfr.integrators.std',
'pyfr.plugins',
'pyfr.quadrules',
'pyfr.readers',
--
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