[pyfr] 06/88: Add a uniform error norm into the PI controller.
Ghislain Vaillant
ghisvail-guest at moszumanska.debian.org
Wed Nov 16 12:05:24 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 22176ceaa9ea41d5570e80449714d05cdd329e82
Author: Freddie Witherden <freddie at witherden.org>
Date: Mon Mar 21 15:03:30 2016 -0700
Add a uniform error norm into the PI controller.
---
pyfr/backends/cuda/blasext.py | 7 +++++--
pyfr/backends/opencl/blasext.py | 7 +++++--
pyfr/backends/openmp/blasext.py | 4 ++--
pyfr/backends/openmp/kernels/errest.mako | 14 ++++++++++----
pyfr/integrators/std/controllers.py | 31 +++++++++++++++++++++++--------
5 files changed, 45 insertions(+), 18 deletions(-)
diff --git a/pyfr/backends/cuda/blasext.py b/pyfr/backends/cuda/blasext.py
index 1d89cf1..30d64dd 100644
--- a/pyfr/backends/cuda/blasext.py
+++ b/pyfr/backends/cuda/blasext.py
@@ -53,7 +53,7 @@ class CUDABlasExtKernels(CUDAKernelProvider):
return CopyKernel()
- def errest(self, x, y, z):
+ def errest(self, x, y, z, *, norm):
if x.traits != y.traits != z.traits:
raise ValueError('Incompatible matrix types')
@@ -62,9 +62,12 @@ class CUDABlasExtKernels(CUDAKernelProvider):
yarr = GPUArray(y.leaddim*y.nrow, y.dtype, gpudata=y)
zarr = GPUArray(z.leaddim*z.nrow, z.dtype, gpudata=z)
+ # Norm type
+ reduce_expr = 'a + b' if norm == 'l2' else 'max(a, b)'
+
# Build the reduction kernel
rkern = ReductionKernel(
- x.dtype, neutral='0', reduce_expr='a + b',
+ x.dtype, neutral='0', reduce_expr=reduce_expr,
map_expr='pow(x[i]/(atol + rtol*max(fabs(y[i]), fabs(z[i]))), 2)',
arguments='{0}* x, {0}* y, {0}* z, {0} atol, {0} rtol'
.format(npdtype_to_ctype(x.dtype))
diff --git a/pyfr/backends/opencl/blasext.py b/pyfr/backends/opencl/blasext.py
index 98faf5d..8fa5a11 100644
--- a/pyfr/backends/opencl/blasext.py
+++ b/pyfr/backends/opencl/blasext.py
@@ -46,16 +46,19 @@ class OpenCLBlasExtKernels(OpenCLKernelProvider):
return CopyKernel()
- def errest(self, x, y, z):
+ def errest(self, x, y, z, *, norm):
if x.traits != y.traits != z.traits:
raise ValueError('Incompatible matrix types')
cnt = x.leaddim*x.nrow
dtype = x.dtype
+ # Norm type
+ reduce_expr = 'a + b' if norm == 'l2' else 'max(a, b)'
+
# Build the reduction kernel
rkern = ReductionKernel(
- self.backend.ctx, dtype, neutral='0', reduce_expr='a + b',
+ self.backend.ctx, dtype, neutral='0', reduce_expr=reduce_expr,
map_expr='pow(x[i]/(atol + rtol*max(fabs(y[i]), fabs(z[i]))), 2)',
arguments='__global {0}* x, __global {0}* y, __global {0}* z, '
'{0} atol, {0} rtol'.format(npdtype_to_ctype(dtype))
diff --git a/pyfr/backends/openmp/blasext.py b/pyfr/backends/openmp/blasext.py
index c97bf4c..8343004 100644
--- a/pyfr/backends/openmp/blasext.py
+++ b/pyfr/backends/openmp/blasext.py
@@ -41,7 +41,7 @@ class OpenMPBlasExtKernels(OpenMPKernelProvider):
return CopyKernel()
- def errest(self, x, y, z):
+ def errest(self, x, y, z, *, norm):
if x.traits != y.traits != z.traits:
raise ValueError('Incompatible matrix types')
@@ -49,7 +49,7 @@ class OpenMPBlasExtKernels(OpenMPKernelProvider):
dtype = x.dtype
# Render the reduction kernel template
- src = self.backend.lookup.get_template('errest').render()
+ src = self.backend.lookup.get_template('errest').render(norm=norm)
# Build
rkern = self._build_kernel(
diff --git a/pyfr/backends/openmp/kernels/errest.mako b/pyfr/backends/openmp/kernels/errest.mako
index 1afdb04..f280904 100644
--- a/pyfr/backends/openmp/kernels/errest.mako
+++ b/pyfr/backends/openmp/kernels/errest.mako
@@ -6,11 +6,17 @@ fpdtype_t
errest(int n, fpdtype_t *__restrict__ x, fpdtype_t *__restrict__ y,
fpdtype_t *__restrict__ z, fpdtype_t atol, fpdtype_t rtol)
{
- fpdtype_t sum = 0.0;
+ fpdtype_t out = 0.0;
- #pragma omp parallel for reduction(+:sum)
+% if norm == 'l2':
+ #pragma omp parallel for reduction(+:out)
for (int i = 0; i < n; i++)
- sum += pow(x[i]/(atol + rtol*max(fabs(y[i]), fabs(z[i]))), 2);
+ out += pow(x[i]/(atol + rtol*max(fabs(y[i]), fabs(z[i]))), 2);
+% else:
+ #pragma omp parallel for reduction(max:out)
+ for (int i = 0; i < n; i++)
+ out = max(out, pow(x[i]/(atol + rtol*max(fabs(y[i]), fabs(z[i]))), 2));
+% endif
- return sum;
+ return out;
}
diff --git a/pyfr/integrators/std/controllers.py b/pyfr/integrators/std/controllers.py
index cefc7da..8875edf 100644
--- a/pyfr/integrators/std/controllers.py
+++ b/pyfr/integrators/std/controllers.py
@@ -87,6 +87,11 @@ class StdPIController(BaseStdController):
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.7)
self._beta = self.cfg.getfloat(sect, 'pi-beta', 0.4)
@@ -105,23 +110,33 @@ class StdPIController(BaseStdController):
@memoize
def _get_errest_kerns(self):
- return self._get_kernels('errest', nargs=3)
+ return self._get_kernels('errest', nargs=3, norm=self._norm)
def _errest(self, x, y, z):
comm, rank, root = get_comm_rank_root()
errest = self._get_errest_kerns()
- # Obtain an estimate for the error
+ # Obtain an estimate for the squared error
self._prepare_reg_banks(x, y, z)
self._queue % errest(self._atol, self._rtol)
- # Reduce locally (element types) and globally (MPI ranks)
- rl = sum(errest.retval)
- rg = comm.allreduce(rl, op=get_mpi('sum'))
-
- # Normalise
- err = math.sqrt(rg / self._gndofs)
+ # L2 norm
+ if self._norm == 'l2':
+ # Reduce locally (element types) and globally (MPI ranks)
+ rl = sum(errest.retval)
+ rg = comm.allreduce(rl, op=get_mpi('sum'))
+
+ # Normalise
+ err = math.sqrt(rg / self._gndofs)
+ # Uniform norm
+ else:
+ # Reduce locally (element types) and globally (MPI ranks)
+ rl = max(errest.retval)
+ rg = comm.allreduce(rl, op=get_mpi('max'))
+
+ # Normalise
+ err = math.sqrt(rg)
return err if not math.isnan(err) else 100
--
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