Merge pull request #684 from evalf/cylinderflow

Cylinderflow

examples/cylinderflow.py

@@ -7,13 +7,68 @@
 # exponentially with radius such that the artificial exterior boundary is
 # placed at large (configurable) distance.
 
-from nutils import mesh, function, solver, util, export, cli, testing
+from nutils import mesh, function, solver, util, export, cli, testing, numeric
 from nutils.expression_v2 import Namespace
 import numpy
 import treelog
 
 
-def main(nelems: int, degree: int, reynolds: float, rotation: float, timestep: float, maxradius: float, seed: int, endtime: float):
+class PostProcessor:
+
+    def __init__(self, topo, ns, timestep, region=4., aspect=16/9, figscale=7.2, spacing=.05):
+        self.ns = ns
+        self.figsize = aspect * figscale, figscale
+        self.bbox = numpy.array([[-.5, aspect-.5], [-.5, .5]]) * region
+        self.bezier = topo.select(function.min(*(ns.x-self.bbox[:,0])*(self.bbox[:,1]-ns.x))).sample('bezier', 5)
+        self.spacing = spacing
+        self.timestep = timestep
+        self.t = 0.
+        self.initialize_xgrd()
+        self.topo = topo
+
+    def initialize_xgrd(self):
+        self.orig = numeric.floor(self.bbox[:,0] / (2*self.spacing)) * 2 - 1
+        nx, ny = numeric.ceil(self.bbox[:,1] / (2*self.spacing)) * 2 + 2 - self.orig
+        self.vacant = numpy.hypot(
+            self.orig[0] + numpy.arange(nx)[:,numpy.newaxis],
+            self.orig[1] + numpy.arange(ny)) > self.ns.R.eval() / self.spacing
+        self.xgrd = (numpy.stack(self.vacant[1::2,1::2].nonzero(), axis=1) * 2 + self.orig + 1) * self.spacing
+
+    def regularize_xgrd(self):
+        # use grid rounding to detect and remove oldest points that have close
+        # neighbours and introduce new points into vacant spots
+        keep = numpy.zeros(len(self.xgrd), dtype=bool)
+        vacant = self.vacant.copy()
+        for i, ind in enumerate(numeric.round(self.xgrd / self.spacing) - self.orig): # points are ordered young to old
+            if all(ind >= 0) and all(ind < vacant.shape) and vacant[tuple(ind)]:
+                vacant[tuple(ind)] = False
+                keep[i] = True
+        roll = numpy.arange(vacant.ndim)-1
+        for _ in roll: # coarsen all dimensions using 3-point window
+            vacant = numeric.overlapping(vacant.transpose(roll), axis=0, n=3)[::2].all(1)
+        newpoints = numpy.stack(vacant.nonzero(), axis=1) * 2 + self.orig + 1
+        self.xgrd = numpy.concatenate([newpoints * self.spacing, self.xgrd[keep]], axis=0)
+
+    def __call__(self, args):
+        x, p = self.bezier.eval(['x_i', 'p'] @ self.ns, **args)
+        logr = numpy.log(numpy.hypot(*self.xgrd.T) / self.ns.R.eval())
+        φ = numpy.arctan2(self.xgrd[:,1], self.xgrd[:,0]) + numpy.pi
+        ugrd = self.topo.locate(self.ns.grid, numpy.stack([logr, φ], axis=1), eps=1, tol=1e-5).eval(self.ns.u, **args)
+        with export.mplfigure('flow.png', figsize=self.figsize) as fig:
+            ax = fig.add_axes([0, 0, 1, 1], yticks=[], xticks=[], frame_on=False, xlim=self.bbox[0], ylim=self.bbox[1])
+            im = ax.tripcolor(*x.T, self.bezier.tri, p, shading='gouraud', cmap='jet')
+            export.plotlines_(ax, x.T, self.bezier.hull, colors='k', linewidths=.1, alpha=.5)
+            ax.quiver(*self.xgrd.T, *ugrd.T, angles='xy', width=1e-3, headwidth=3e3, headlength=5e3, headaxislength=2e3, zorder=9, alpha=.5, pivot='tip')
+            ax.plot(0, 0, 'k', marker=(3, 2, self.t*self.ns.rotation.eval()*180/numpy.pi-90), markersize=20)
+            cax = fig.add_axes([0.9, 0.1, 0.01, 0.8])
+            cax.tick_params(labelsize='large')
+            fig.colorbar(im, cax=cax)
+        self.t += self.timestep
+        self.xgrd += ugrd * self.timestep
+        self.regularize_xgrd()
+
+
+def main(nelems: int, degree: int, reynolds: float, rotation: float, radius: float, timestep: float, maxradius: float, seed: int, endtime: float):
     '''
     Flow around a cylinder.
 
@@ -28,6 +83,8 @@ def main(nelems: int, degree: int, reynolds: float, rotation: float, timestep: f
          less.
        reynolds [1000]
          Reynolds number, taking the cylinder radius as characteristic length.
+       radius [.5]
+         Cylinder radius.
        rotation [0]
          Cylinder rotation speed.
        timestep [.04]
@@ -51,10 +108,12 @@ def main(nelems: int, degree: int, reynolds: float, rotation: float, timestep: f
     ns.δ = function.eye(domain.ndims)
     ns.Σ = function.ones([domain.ndims])
     ns.uinf = function.Array.cast([1, 0])
-    ns.r = .5 * function.exp(elemangle * geom[0])
     ns.Re = reynolds
-    ns.phi = geom[1] * elemangle  # add small angle to break element symmetry
-    ns.x_i = 'r (cos(phi) δ_i0 + sin(phi) δ_i1)'
+    ns.grid = geom * elemangle
+    ns.R = radius
+    ns.r = 'R exp(grid_0)'
+    ns.φ = 'grid_1'
+    ns.x_i = 'r (cos(φ) δ_i0 + sin(φ) δ_i1)'
     ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
     J = ns.x.grad(geom)
     detJ = function.determinant(J)
@@ -66,7 +125,7 @@ def main(nelems: int, degree: int, reynolds: float, rotation: float, timestep: f
     ns.N = 10 * degree / elemangle  # Nitsche constant based on element size = elemangle/2
     ns.nitsche_i = '(N v_i - (∇_j(v_i) + ∇_i(v_j)) n_j) / Re'
     ns.rotation = rotation
-    ns.uwall_i = '0.5 rotation (-sin(phi) δ_i0 + cos(phi) δ_i1)'
+    ns.uwall_i = '0.5 rotation (-sin(φ) δ_i0 + cos(φ) δ_i1)'
 
     inflow = domain.boundary['outer'].select(-ns.uinf.dotnorm(ns.x), ischeme='gauss1')  # upstream half of the exterior boundary
     sqr = inflow.integral('Σ_i (u_i - uinf_i)^2' @ ns, degree=degree*2)
@@ -82,35 +141,15 @@ def main(nelems: int, degree: int, reynolds: float, rotation: float, timestep: f
     res += domain.integral('q ∇_k(u_k) dV' @ ns, degree=9)
     uinertia = domain.integral('v_i u_i dV' @ ns, degree=9)
 
-    bbox = numpy.array([[-2, 46/9], [-2, 2]])  # bounding box for figure based on 16x9 aspect ratio
-    bezier0 = domain.sample('bezier', 5)
-    bezier = bezier0.subset((bezier0.eval((ns.x-bbox[:, 0]) * (bbox[:, 1]-ns.x)) > 0).all(axis=1))
-    interpolate = util.tri_interpolator(bezier.tri, bezier.eval(ns.x), mergetol=1e-5)  # interpolator for quivers
-    spacing = .05  # initial quiver spacing
-    xgrd = util.regularize(bbox, spacing)
+    postprocess = PostProcessor(domain, ns, timestep)
 
     with treelog.iter.plain('timestep', solver.impliciteuler('u:v,p:q', residual=res, inertia=uinertia, arguments=args0, timestep=timestep, constrain=cons, newtontol=1e-10)) as steps:
         for istep, args in enumerate(steps):
 
-            t = istep * timestep
-            x, u, normu, p = bezier.eval(['x_i', 'u_i', 'sqrt(u_i u_i)', 'p'] @ ns, **args)
-            ugrd = interpolate[xgrd](u)
-
-            with export.mplfigure('flow.png', figsize=(12.8, 7.2)) as fig:
-                ax = fig.add_axes([0, 0, 1, 1], yticks=[], xticks=[], frame_on=False, xlim=bbox[0], ylim=bbox[1])
-                im = ax.tripcolor(x[:, 0], x[:, 1], bezier.tri, p, shading='gouraud', cmap='jet')
-                import matplotlib.collections
-                ax.add_collection(matplotlib.collections.LineCollection(x[bezier.hull], colors='k', linewidths=.1, alpha=.5))
-                ax.quiver(xgrd[:, 0], xgrd[:, 1], ugrd[:, 0], ugrd[:, 1], angles='xy', width=1e-3, headwidth=3e3, headlength=5e3, headaxislength=2e3, zorder=9, alpha=.5)
-                ax.plot(0, 0, 'k', marker=(3, 2, t*rotation*180/numpy.pi-90), markersize=20)
-                cax = fig.add_axes([0.9, 0.1, 0.01, 0.8])
-                cax.tick_params(labelsize='large')
-                fig.colorbar(im, cax=cax)
-
-            if t >= endtime:
-                break
+            postprocess(args)
 
-            xgrd = util.regularize(bbox, spacing, xgrd + ugrd * timestep)
+            if istep * timestep >= endtime:
+                break
 
     return args0, args
 
@@ -131,7 +170,7 @@ def main(nelems: int, degree: int, reynolds: float, rotation: float, timestep: f
 class test(testing.TestCase):
 
     def test_rot0(self):
-        args0, args = main(nelems=6, degree=3, reynolds=100, rotation=0, timestep=.1, maxradius=25, seed=0, endtime=.05)
+        args0, args = main(nelems=6, degree=3, reynolds=100, radius=.5, rotation=0, timestep=.1, maxradius=25, seed=0, endtime=.05)
         with self.subTest('initial condition'):
             self.assertAlmostEqual64(args0['u'], '''
                 eNoNzLEJwkAYQOHXuYZgJyZHLjmMCIJY2WcEwQFcwcYFbBxCEaytzSVnuM7CQkEDqdVCG//uNe/rqEcI
@@ -148,7 +187,7 @@ def test_rot0(self):
                 wL9FwkabQsJJTbc2ubJHw7ZvRq/qITA=''')
 
     def test_rot1(self):
-        args0, args = main(nelems=6, degree=3, reynolds=100, rotation=1, timestep=.1, maxradius=25, seed=0, endtime=.05)
+        args0, args = main(nelems=6, degree=3, reynolds=100, radius=.5, rotation=1, timestep=.1, maxradius=25, seed=0, endtime=.05)
         with self.subTest('initial condition'):
             self.assertAlmostEqual64(args0['u'], '''
                 eNoNzLEJwkAYQOHXuYZgJyZHLjmMCIJY2WcEwQFcwcYFbBxCEaytzSVnuM7CQkEDqdVCG//uNe/rqEcI

nutils/topology.py

@@ -1787,6 +1787,9 @@ def refined(self):
         groups = [{name: topo.refined if isinstance(topo, TransformChainsTopology) else topo for name, topo in groups.items()} for groups in (self.vgroups, self.bgroups, self.igroups, self.pgroups)]
         return self.basetopo.refined.withgroups(*groups)
 
+    def locate(self, geom, coords, **kwargs):
+        return self.basetopo.locate(geom, coords, **kwargs)
+
 
 class OppositeTopology(TransformChainsTopology):
     'opposite topology'
@@ -1840,7 +1843,7 @@ def StructuredLine(space, root: transform.stricttransformitem, i: types.strictin
 class StructuredTopology(TransformChainsTopology):
     'structured topology'
 
-    __slots__ = 'root', 'axes', 'nrefine', 'shape', '_bnames'
+    __slots__ = 'root', 'axes', 'nrefine', 'shape', '_bnames', '_asaffine_geom', '_asaffine_retval'
     __cache__ = 'connectivity', 'boundary', 'interfaces'
 
     @types.apply_annotations
@@ -2237,37 +2240,42 @@ def refined(self):
         axes = [axis.refined for axis in self.axes]
         return StructuredTopology(self.space, self.root, axes, self.nrefine+1, bnames=self._bnames)
 
-    def locate(self, geom, coords, *, tol, eps=0, weights=None, skip_missing=False, arguments=None, **kwargs):
+    def locate(self, geom, coords, *, tol=0, eps=0, weights=None, skip_missing=False, arguments=None, **kwargs):
         coords = numpy.asarray(coords, dtype=float)
         if geom.ndim == 0:
             geom = geom[_]
             coords = coords[..., _]
         if not geom.shape == coords.shape[1:] == (self.ndims,):
             raise Exception('invalid geometry or point shape for {}D topology'.format(self.ndims))
-        arguments = dict(arguments or ())
-        geom0, scale, e = self._asaffine(geom, arguments)
-        if e > tol:
-            return super().locate(geom, coords, eps=eps, tol=tol, weights=weights, skip_missing=skip_missing, arguments=arguments, **kwargs)
-        log.info('locate detected linear geometry: x = {} + {} xi ~{:+.1e}'.format(geom0, scale, e))
-        return self._locate(geom0, scale, coords, eps=eps, weights=weights, skip_missing=skip_missing)
+        if tol or eps:
+            if arguments:
+                geom0, scale, error = self._asaffine(geom, arguments)
+            elif geom is getattr(self, '_asaffine_geom', None):
+                log.debug('locate found previously computed affine values')
+                geom0, scale, error = self._asaffine_retval
+            else:
+                self._asaffine_geom = geom
+                geom0, scale, error = self._asaffine_retval = self._asaffine(geom, {})
+            if all(error <= numpy.maximum(tol, eps * scale)):
+                log.debug('locate detected linear geometry: x = {} + {} xi ~{}'.format(geom0, scale, error))
+                return self._locate(geom0, scale, coords, eps=eps, weights=weights, skip_missing=skip_missing)
+        return super().locate(geom, coords, eps=eps, tol=tol, weights=weights, skip_missing=skip_missing, arguments=arguments, **kwargs)
 
     def _asaffine(self, geom, arguments):
         # determine geom0, scale, error such that geom ~= geom0 + index * scale + error
-        funcsp = self.basis('std', degree=1, periodic=())
-        verts = numeric.meshgrid(*map(numpy.arange, numpy.array(self.shape)+1)).reshape(self.ndims, -1)
-        index = (funcsp * verts).sum(-1)
-        # strategy: fit an affine plane through the minima and maxima of a
-        # uniform sample, and evaluate the error as the largest difference on
-        # the remaining sample points
-        geom_, index_ = self.sample('uniform', 2 + (1 in self.shape)).eval((geom, index), **arguments)
-        imin = geom_.argmin(axis=0)
-        imax = geom_.argmax(axis=0)
-        R = numpy.arange(self.ndims)
-        scale = (geom_[imax,R] - geom_[imin,R]) / (index_[imax,R] - index_[imin,R])
-        geom0 = geom_[imin,R] - index_[imin,R] * scale # geom_[im..,R] = index_[im..,R] * scale + geom0
-        error = numpy.linalg.norm(geom0 + index_ * scale - geom_, axis=1).max()
-
-        return geom0, scale, error
+        n = 2 + (1 in self.shape) # number of sample points required to establish nonlinearity
+        sampleshape = numpy.multiply(self.shape, n) # shape of uniform sample
+        geom_ = self.sample('uniform', n).eval(geom, **arguments) \
+            .reshape(*self.shape, *[n] * self.ndims, self.ndims) \
+            .transpose(*(i+j for i in range(self.ndims) for j in (0, self.ndims)), self.ndims*2) \
+            .reshape(*sampleshape, self.ndims)
+        # strategy: fit an affine plane through the minima and maxima of a uniform sample,
+        # and evaluate the error as the largest difference on the remaining sample points
+        xmin, xmax = geom_.reshape(-1, self.ndims)[[0, -1]]
+        dx = (xmax - xmin) / (sampleshape-1) # x = x0 + dx * (i + .5) => xmax - xmin = dx * (sampleshape-1)
+        for idim in range(self.ndims):
+            geom_[...,idim] -= xmin[idim] + dx[idim] * numpy.arange(sampleshape[idim]).reshape([-1 if i == idim else 1 for i in range(self.ndims)])
+        return xmin - dx/2, dx * n, numpy.abs(geom_).reshape(-1, self.ndims).max(axis=0)
 
     def _locate(self, geom0, scale, coords, *, eps=0, weights=None, skip_missing=False):
         mincoords, maxcoords = numpy.sort([geom0, geom0 + scale * self.shape], axis=0)

tests/test_topology.py

@@ -959,7 +959,7 @@ def test_missing_argument(self):
     @parametrize.enable_if(lambda etype, mode, **kwargs: etype == 'square' and mode != 'nonlinear')
     def test_detect_linear(self):
         target = numpy.array([(.2, .3)])
-        with self.assertLogs('nutils', level='INFO') as cm:
+        with self.assertLogs('nutils', level='DEBUG') as cm:
             self.domain.locate(self.geom, target, eps=1e-15, tol=1e-12, arguments=dict(scale=.123))
         self.assertRegex(cm.output[0], 'locate detected linear geometry')